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INTRODUCTION 



Semiconductor devices operate by controlling the flow of electrons and holes through a device, and our understanding 
of charge carrier transport has both benefited from and driven t he deve lopment of semiconductor devices. When 
' William Shockley wrote " Electrons and Holes in Semiconductors" [Sho50l | , semiconductor physics was at the frontier 
[ of research in condensed matter physics. Over the years, the essential concepts were clarified and simplified into 
• the working knowledge of device engineers. The treatment of electrons and holes as semiclassical particles with an 
1) ' effective mass was usually adequate. Electronic devices were made of materials (e.g. silicon, gallium arsenide, etc.) 
with properties (e.g. bandgap, effective mass, etc.) that could be looked up. For most devices, the engineer's drift- 
diffusion equation provided a simple, but adequate description of carrier transport. Today things are changing. Device 
dimensions have shrunk to the nanoscale. The properties of materials can be engineered by intentional strain and size 
effects due to quantum confinement. Devices contain a countable number of dopants and are sensitive to structure at 
the atomistic scale. In addition to familiar devices like the MOSFET, which have been scaled to nanometer dimensions, 
new devices built from carbon nanotubes, semiconductor nanowires, and organic molecules are being explored. Device 
engineers will need to learn to think about devices differently. To describe carrier transport in nanoscale devices, 
engineers must learn how to think about charge carriers as quantum mechanical entities rather than as semiclassical 
particles, and they must learn how to think at the atomistic scale rather than at a continuum one. Our purpose in 
this chapter is to provide engineers with an introduction to the non-equilibrium Green's function (NEGF) approach, 
which provides a powerful conceptual tool and a practical analysis method to treat small electronic devices quantum 
mechanically and atomistically. We first review the basis for the traditional, semiclassical description of carriers that 
has served device engineers for more than 50 years in section [XXJ We then describe why this traditional approach 
\ loses validity at the nanoscale. Next, we describe semiclassical ballistic transport in section Hill and the Landauer- 
f^i . Buttikcr approach to phase coherent quantum transport in section IIVI Realistic devices include interactions that 
t-H ' break quantum mechanical phase and also cause energy relaxation. As a result, transport in nanodevices are between 
^^O , diffusive and phase coherent. We introduce the non equilbrium Green's function (NEGF) approach, which can be used 
to model devices all the way from ballistic to diffusive limits in section [V] This is followed by a summary of equations 
that are used to model a large class of layered structures such as nanotransistors, carbon nanotubes and nanowires 
in section IVI1 An application of the NEGF method in the ballistic and scattering limits to silicon nanotransistors is 
discussed in sections IVHI and IVIIII respectively We conclude with a summary in section IIXI The Dyson's equations 
. and algorithms to solve for the Green's functions of layered structures are presented in appendices B and C. These 
appendices can be left out by a reader whose aim is to gain a basic understanding of the NEGF approach to device 
O ■ modeling. 
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II. SEMICLASSICAL TRANSPORT: DIFFUSIVE 



Electrical engineers have commonly treated electrons as semiclassical particles that move through a device under 
the influence of an electric field and random scattering potentials. As sketched in Fig. 1, electrons move along a 
trajectory in phase space (position and momentum space). In momentum space, the equation of motion looks like 
Newton's Law for a classical particle 

dfrk „ „ „ 

— = -V r E c (r,t) , (1) 
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FIG. 1: Carrier trajectories in phase space showing free flights along a trajectory interrupted by scattering events that begin 
another free flight. 



where E(k) describes the bandstructure of the semiconductor. The right hand side of eqn. @ is simply the velocity 
of the semiclassical particle, and in the simplest case it is just hk/m*. By solving eqns. ([1]) and ([2]), we trace the 
trajectory of a carrier in phase space as shown in Fig. [TJ 

Equations (Q]) and @ describe the ballistic transport of semiclassical carriers. In practice, carriers frequently scatter 
from various perturbing potentials (defects, ionized impurities, lattice vibrations, etc.). The result is that carriers hop 
from one trajectory in phase space to another as shown in Fig. 1. The average distance between scattering events, the 
mean- free-path, I, has (until recently) been much smaller than the critical dimensions of a device. Carriers undergo 
a random walk through a device with a small bias in one direction imposed by the electric field. To describe this 
scattering dominated (so-called diffusive) transport, we should add a random force (Fs(r,t)) to the right hand side 
of eqn. (1) 



It is relatively easy to solve eqns. (fTJ) and ((3J numerically. One solves the equations of motion, eqns. (1) and (2) to 
move a particle through phase space. Random numbers are chosen to mimic the scattering process and occasionally 
kick a carrier to another trajectory. By averaging the results for a large number of simulated trajectories, these so- 
called Monte Carlo techniques provide a rigorous, though computationally demanding description of carrier transport 
in devices as described in [Jac89]. 

Device engineers are primarily interested in average quantities such as the average electron density, current density, 
etc. (There are some exceptions; noise is important too.) Instead of simulating a large number of particles, we can 
ask: What is the probability that a state at position, r, with momentum hk, is occupied at time i? The answer is 
given by the distribution function, f(r,k,t), which can be computed by averaging the results of a large number of 
simulated trajectories. Alternatively, we can adopt a collective viewpoint instead of the individual particle viewpoint 
and formulate an equation for /(r, k,t). The result is known as the Boltzmann Transport Equation (BTE), 



where k is the crystal momentum and Eq is the bottom of the conduction band. In position space, the equation of 
motion is 




dhk 



V r Ec(r 7 t) + F s (f 7 t) , 



(3) 



(4) 



where E is the electric field, and Cf describes the effects of scattering. In equilibrium f(r,k,t) is simply the Fermi 
function, but in general, we need to solve eqn. ([?]) to find f. Once /(f, k, t) is known, quantities of interest to the device 
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engineer are readily found. For example, to find the average electron density in a volume CI centered at position, r , 
we simply add up the probability that all of the states in O are occupied and divide by the volume, 



(5) 



Similarly, we find, 



Current Density J(r,t) = i ^~^(— o)vf(r, k,t) , 

k 

Kinetic Energy Density W(r, t) = ^ E(k)f(f, k, t) , 

k 

Energy Current Density Ji;(r, t) = — ^ E(k)vf(r, k, t) 



(6) 
(7) 
(8) 



where q > is the absolute value of the electron charge. This approach provides a clear and fairly rigorous description 
of semiclassical carrier transport, but solving the six dimensional BTE is enormously difficult. One might ask if we 
can't just find a way to solve directly for the quantities of interest in eqns. (JSJ) - ©. The answer is yes, but some 
simplifying assumptions are necessary. 

Device engineers commonly describe carrier transport by few low order moments of the Boltzmann transport equa- 
tion, eqn. A mathematical prescription for generating mo ment equ ations exists, but to formulate them in a 
tractable manner, numerous simplifying assumptions are required EunOfJ. Moment equations provide a phenomeno- 
logical description of transport that gives insight and quantitative results when properly calibrated. 

The equation for the zeroth moment of f(f, k,t) gives the well known continuity equation for the electron density, 
n(r,t), 



dn(f, t) 
dt 



where F n is the electron flux, 



(9) 



(10) 



G n the electron generation rate, and R n the electron recombination rate. Equation ^ states that the electron density 
at a location increases with time if there is a net flux of electrons into the region (as described by the first term, 
minus the divergence of the electron flux) or if carriers are being generated there. Recombination causes the electron 
density to decrease with time. Any physical quantity must obey a conservation law like eqn. ([9]). 

The equation for the first moment of f(r, k,t) gives the equation for average current density (eqn. ([6])) projection 
on the a;-axis, 



d I 

UtJ nj 



2q dW x , 
m dx 



nq 



E x 
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(11) 



Each term on the right hand side of eqn. (fTT|) is analogous to the corresponding terms in eqn. ([9|). The current 
typically changes slowly on the scale of the momentum relaxation time, r m , (typically a sub-picosccond time) so the 
time derivative can be ignored and eqn. pip solved for 



T TP I 2 dW *- 

nqfini^x -t- —fi n — - — 

6 ax 



where 



m* 



(12) 



(13) 



is the electron mobility, and we have assumed equipartition of energy so that W xx = W/3, where W is the total 
kinetic energy density. This assumption can be justified when there is a lot of isotropic scattering, which randomizes 
the carrier velocity. Eqn. (Ti"2")) is a drift-diffusion equation; it says that electrons drift in electric fields and diffuse 
down kinetic energy gradients. Near equilibrium, 



W = -nk B T , 



(14) 
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FIG. 2: The average velocity vs. electric field for electrons in bulk silicon at room temperature. 



so when T is uniform, cqn. (j 1 2[) becomes 



dfi dn 
Jnx - nqn n E x + k B T^ n — = nq/j, n E x + 1 D n-r-, (15) 

the drift-diffusion equation. By inserting eqn. (fT5)) in the electron continuity equation, cqn. ([9]), we get an equation 
for the electron density that can be solved for the electron density within a devi ce. This is the traditional and still 
most common approach for describing transport in semiconductor devices |Pie87l | . 

Since most devices contain regions with high electric fields, the assumption that W = inksT/2 is not usually a 
good one. The carrier energy enters directly into the second term of the transport equation, eqn. (fT2|) . but it also 
enters indirectly because the mobility is energy dependent. To treat transport more rigoursly, we need an equation 
for the electron energy. 

The second moment of f(r, k, t) gives the carrier energy den sity, W (r, t) according to eqn ([5]). The second moment 
of the BTE gives a continuity equation for the energy density LunOOj ] 



dw dj w w-w 

-xr = — j J nxE x , (16) 

at dx te 

where Wo is the equilibrium energy density and te the energy relaxation time. Note that the energy relaxation time 
is generally longer than the momentum relaxation time because phonon energies are small so that it takes several 
scattering events to thermalize an energetic carrier but only one to randomize its momentum. 

To solve eqn. (|I6j) . we need to specify the energy current. The third moment of f(r,k,t) gives the carrier energy 
flux, Jwir, t) according to eqn. ||SJ). The third moment of the BTE gives a continuity equation for the energy flux, 

Jw = W, E E X + ^fD. , (17) 
dx 

where he an( l De arc appropriate energy transport mobility and diffusion coefficient [LunO0| | . 

Equations (J9j> , (fl5|) . (|16|) and (fT7|) can now be solved seli-consistently to simulate carrier transport. Fig. ([2]) sketches 
the result for uniformly doped, bulk silicon with a constant electric field. At low electric fields, W ~ 3nfcsT/2, and 
we find that < v x >~ [i n E x . For electric fields above ~ 10 V/ cm, the kinetic energy increases, which increases 
the rate of scattering and lowers the mobility so that at high fields the velcoity saturates at ~ 10 7 cm/ s. In a bulk 
semiconductor, there is a one-to-one relation between the magnitude of the electric field and the kinetic energy, so 
the mobility and diffusion coefficient can be parametrized as known functions of the local electric field. The result is 
that for bulk semiconductors or for large devices in which the electric field changes slowly, there is no need to solve all 
four equations; we need to solve the carrier continuity and drift-diffusion equations with ficld-dcpcndcnt parameters. 

Electric fields above lQ A V/cm arc common in nanoscale devices. This is certainly high enough to cause velocity 
saturation in the bulk, but in a short, high field region, transients occur. Fig. [3] illustrates what happens for a 
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FIG. 3: The average steady-state velocity and kinetic energy vs. postion for electrons injected into a short slab of silicon with 
low-high-low electric field profile. 



hypothetical situation in which the electric field abrubtly jumps from a low value to a high value and then back to 
a low value again. Electrons injected from the low field region are accelerated by the high electric field, but energy 
relaxation times are longer than momentum relaxation times, so the eneregy is slow to respond. The result is that 
the mobility is initially high (even though the electric field is high), so the velocity can be higher than the saturated 
value shown in Fig. O As the kinetic energy increases, however, scattering increases, the mobility drops, and the 
velocity eventually decreases to ~ 10 cm/ s, the saturated velocity for electrons in bulk silicon. The spatial width 
of the transient is roughly 100 nm; modern devices frequently have dimensions on this order, and strong velocity 
overshoot should be expected. 

The example shown in Fig. [3] demonstrates that it is better to think of the mobility and diffusion coefficient as 
functions of the local kinetic energy rather than the local electric field. What this means is that eqns. (O, (JTSJ) , (fTB)) and 
(fT7|) should all be solved self-consistently to simulate carrier transport in small devices. Device simulation programs 
commonly permit two options: 1) the solution of eqns. ([9]) and (|15[) self-consistently with the Poisson equation using 
mobilities and diffusion coefficients that depend on the local electric field (the so-called drift-diffusion approach) or 2) 
the solution of eqns. (fTS)) . (fTrj|) and (Q~7|) self-consistently with the Poisson equation using mobilities and diffusion 
coefficients that depend on the local kinetic energy (so-called energy transport or hydrodynamic approaches). Solving 
the four equations self-consistently is more of a computational burden, but it is necessary for when the electric field 
changes rapidly on the scale of a mean-free-path for scattering. Actually, numerous sim plifying assumptions are also 



necessary to even write the current and energy flux equations as eqns. (|15p and (| 1 7[) LunOOl ] . The most rigorous 
(and computationally demanding) simulations (so-called Monte Carlo simulations) go back to the individual particle 
picture and track carriers trajectories according to eqns. ([2]) and ([3]). Drift-diffusion and energy transport approaches 
for treating carrier transport in semiconductor devices have two things in common; the first is the assumption that 
carriers can be treated as semiclassical particles and the second is the assumption that there is a lot of scattering. 
Both of these assumptions arc losing validity as devices shrink. 



III. SEMICLASSICAL TRANSPORT: BALLISTIC 



Consider the "device" sketched in Fig. UK a), which consists of a ballistic region attached to two contacts. The left 
contact (source) injects a thermal equilibrium flux of carriers into the device; some carriers reflect from the potential 
barriers within the device, and the rest transmit across and enter the right contact (drain). A similar statement 



6 



applies to the drain contact. The source and drain contacts arc assumed to be perfect absorbers, which means that 
carriers impinging them from the device travel without reflecting back into the device. To compute the electron 
density, current, average velocity, etc. within the device, we have two choices. The first choice treats the carriers as 
semiclassical particles and the Boltzmann equation is solved to obtain the distribution function f(f, k, t) as discussed 
in the previous section. The second choice treats the carriers quantum mechanically as discussed in the next section. 
In this section, we will use a semiclassical description in which the local density-of-states within the device is just 
that of a bulk semiconductor, but shifted up or down by the local electrostatic potential. This approximation works 
well when the electrostatic potential does not vary too rapidly, so that quantum effects can be ignored. To find how 
the k-states within the ballistic device are occupied, we solve the Boltzmann Transport Equation, eqn. (4). Because 
the device is ballistic, there is no scattering, and Cf = 0. It can be shown Lun0fj| | that the solution to the BTE with 
Cf — is any function of the electron's total energy, 

E = E c (x)+E(k) (18) 

where, Ec(x) is the conduction band minimum versus position and E(k) is the band structure for the conduction 
band. We know that under equilibrium conditions sketched in Fig. 4(b), the proper function of total energy is the 
Fermi function, 

f( E ) = 7- — yra- , (19) 
i + «cp(iy^) 

where the Fermi level, Ep, and temperature, T are constant in equilibrium. 

Now consider the situation in Fig. 4c where a drain bias has been applied to the ballistic device. Although two 
thermal equilibrium fluxes are injected into the device, it is now very far from equilibrium. Since scattering is what 
drives the system to equilibrium, the ballistic device is as far from equilibrium as it can be. Nevertheless, for the 
ballistic device, the relevant steady-state Boltzmann equation is the same equation as in equilibrium. The solution is 
again a function of the carrier's total kinetic energy. At the contacts, we know that the solution is a Fermi function, 
which specifics the functional dependence on energy. For the ballistic device, therefore, the probability that a k-state 
is occupied is given by an equilibrium Fermi function. The only difficulty is that we have two Fermi levels, so we need 
to decide which one to use. 

Return again to Fig. 4c and consider how to fill the states at x = x\. We know that the probability that a k-state 
is occupied is given by a Fermi function, so we only need to decide which Fermi level to use for each k-state. For the 
positive k-states with energy above E^op i the top of the energy barrier, the states can only have been occupied by 
injection from the source, so the appropriate Fermi level to use is the source Fermi level. Similarly, negative k-states 
with energy above Etop can only be occupied by injection from the drain, so the appropriate Fermi level to use 
is the drain Fermi level. Finally, for k-states below EtoPi both positive and negative velocity states are populated 
according to the drain Fermi level. The negative velocity k-states are populated directly by injection from the drain 
and the positive k-states are populated when negative velocity carriers reflect from the potential barrier. 

Ballistic transport can be viewed as a special kind of equilibrium. Each k-state is in equilibrium with the contact 
from which it was populated. Using this reasoning, one can compute the distribution function and any moment 
of it (e.g. carrier density, carrier velocity, etc.) at any location within the device. Figure 5 shows that computed 
distribution function in a ballistic nanoscale MOSFET under high gate and drain bias Rhe02| . A strong ballistic peak 



develops as carriers are injected from the source are accelerated in the high electric field near the drain. Each k-state is 
in equilibrium with one of the two contacts, but the overall carrier distribution is very different from the equilibrium 
Fermi-Dirac distribution. When scattering dominates, carriers quickly lose their "memory" of which contact they 
were injected from, but for ballistic transport there are two separate streams of carriers; one injected from the source 
and one from the drain. 

To evaluate the electron density vs. position within the ballistic device, we should compute a sum like eqn. l[5]). 
but we must do two sums, one for the states filled from the left contact and another for the states filled from the 
right contact, 

n(x) = &L(x)f L (E) + &R(x)f R (E) , (20) 

where /l and fn are the equilibrium Fermi functions of contacts L and R and &l,r(x) is a function that selects out 
the k-states at position, x, that can be filled by contact L or R according to the procedure summarized in Fig. 4. It 
is often convenient to do the integrals in energy space rather than in k-space in which case, eqn. (|20p becomes 



n(x) = J dE LDOS L {x,E)f L (E) + J LDOS R {x,E)f R (E) 



(21) 
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FIG. 4: Sketch of a ballistic device with two contacts that function as reservoirs of thermal equilibrium carriers, (a) Device 
and the two contacts, (b) Energy band diagram under equilibrium conditions (Vn = 0). (c) Energy band diagram under bias 
(Vb > 0). 



where LDOSl.r( x , E) is the local density of states at energy E, fillablc from contact L or R. For diffusive transport, 
we deal with a single density-of-states and fill it according to a since quasi-Fermi level, but for ballistic devices, the 
density of states separates into parts fillable from each contact. 

The current flowing from source to drain (drain to source) contact is simply the transmission probability T(E), 
times the Fermi function of the source (drain) contact. The net current flowing in the device is then, 

I = ~2 J dET(E) [f L (E) - fji(E)] . (22) 

For the semiclassical example of Fig. |H T(E)=0 for E < E T op and T(E)=1 for E > E T op- 

IV. PHASE COHERENT QUANTUM TRANSPORT: THE LANDAUER-BUTTIKER FORMALISM 

Quantum mechanically the electron is a wave and the wave function \l/(r) is obtained by solving Schrodinger's 
equation, 



V n (r) = E n V n (r) , (23) 
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FIG. 5: All wave funcitons in the device can be represented as left incident ), right incident or localized states 

(< C) )- 

where £"„ is the energy. Consider a device connected to two contacts as shown in Fig. where the contacts are 
assumed to have a constant electrostatic potential. In a manner identical to the discussion of semiclassical modeling in 
the previous section, where we considered corpuscular electrons incident from the left and right contacts, in quantum 
mechanical modeling, we need to consider electron waves incident from the left and right contacts. The electron wave 
function in device region D can be thought to arise from: 

• Waves incident from left lead of the form e lkx , which have transmitted and reflected components te lk x and re lkx 
in the right and left leads respectively. The wave function in the device region D due to this wave is represented 
by *£\ 

• Waves incident from right lead of the form e~ lkx , which have transmitted and reflected components te~ lk x and 
re lkx in the right and left leads respectively. The wave function in the device region D due to this wave is 
represented by \E'^f' . 

• States localized in Device region represented by ^>^ C K Localized and quasi-localized states are filled up by 
scattering due to electron-phonon and electron-electron interaction. We will assume here that localized states 
are absent. 

The Landauer-Buttikcr approach expresses the expectation value of an operator in terms of the left and right 
incident electrons from the contacts and their distribution functions. The expectation value of operator Q is: 

Q = E [< ^dU^d > h(E)+ < tfjf^l*^ > fB.(E)] . (24) 

k,s 

Here the summation is performed over the momentum, k, and the spin, s, states. Since in this chapter we do not 
consider any spin-depnednt phenomena, the spin summation trnaslates into a factor of 2. We will distinguish this 
factor by keeping in front of sums and integrals. Eq. p4|) has contributions from two physically different sources. The 



first term corresponds to contribution from waves incident from the left contact ) & t energy E, weighted by the 
Fermi factor of the left contact (/l). The second term corresponds to waves incident from the right contact (^^) 
weighted by the Fermi factor of the right contact (/r). More generally, if region D is connected to a third contact G, 
then the expectation value of operator Q is, 

Q = Y,\< ^IQI*^ > ME)+ < > f R (E)+ < ^ { n } \Q\^D } > Ig(E)] , (25) 
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where corresponds to the wave function in the Device due to waves incident from contact G and fa is the Fermi 
factor of contact G. Using eqn. the contribution to electron (n) and current (J) densities at x in the device 

region D are given by, 



n(x) = (*)\ 2 fL{E) + \^\x)\ 2 f R (E) 

k.s 



and 



eh 
2mi 



c.c 



(26) 



(27) 



The quantum mechanical density of states at energy E due to waves incident from the left (LDOSl) and right 
{LDOSr) are: 



LDOS L (x,E) = 2 l*^)! 3 
k states with energy E incident from the left contact 

LDOS R ( X> E) = 2 J2 \*£H*)\ 2 

k states with energy E incident from the right contact 

Then the electron density can be written in the same form as eqn. (1211) . 

n(x) - 



LDOS L (x,E)f L (E)dE + J LDOS R (x,E)f R (E)dE 



(28) 
(29) 

(30) 



except that the expressions for the local density of states are different. Similarly, eqn. (|2~?| can be expressed in a 
form identical to eqn. (|2"2"|) . The above formalism can be extended to calculate noise (shot and Johnson- Nyquist) 
in nanodevices [Buttikcr92]. Device modeling in the phase coherent limit involves solving Schrodingcr's equation to 
obtain the electron density self-consistcntly with Poisson's equation. 



V. QUANTUM TRANSPORT WITH SCATTERING: THE NEED FOR GREEN'S FUNCTIONS 

The description in the previous section is valid only in the phase coherent limit. The terminology "phase coherent" 
refers to a deterministic evolution of both the amplitude and phase of x & n (r) as given by Schrodinger's equation. 
The quantum mechanical wave function evolves phase coherently only in the presence of rigid scatterers, a common 
example of which is the electrostatic potential felt by an electron in the device. The wave function of an electron 
loses phase coherence due to scatterers which have an internal degree of freedom such as phonons. Phase incoherent 
scattering involves irreversible loss of phase information to phonon degrees of freedom. Naturally, including loss of 
phase information is important when device dimensions become comparable to the scattering lengths due to phonons 
and other phase breaking mechanisms. Accurate modeling of nanodevices should have the ability to capture: 

• Interference effects 

• Quantum mechanical tunneling 

• Discrete energy levels due to confinement in 2D and 3D device geometries 

• Scattering mechanisms (electron-phonon, electron-electron). 

The first three effects can be modeled by solving Schodinger's equation in a rigid potential as discussed in section 
IIVI While in the scmiclassical device modeling, the Boltzmann equation accounts for the energy and momentum 
relaxation due to scattering mechanisms, in quantum mechanical device modeling, the NEGF approach is necessary 
to account for energy, momentum and quantum mechanical phase relaxation. 

In the remainder of this sec tion, we explain the NEGF approach in the phase coherent limit by starting from 
Schrodinger's equation Dat96l |. We will start by an explanation of the tight binding Hamiltonian and relate this 



Hamiltonian to a device with open boundary conditions (section IV Ap . The open boundary conditions lead to an 
infinite dimensional matrix. We will describe a procedure to fold the effect of the open boundaries into the finite 
device region in section [VBI This will allow us to deal with small matrices where the open boundaries are modeled 
by contact self energies. The Green's functions, self energies and their relationship to current and electron density are 
derived in section IV CI Then in section IV Dl we extend the discussion in section IV CI to one include electron-phonon 
interaction, which is where the NEGF approach is really essential. 
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FIG. 6: A one dimensional device connected to two semi-infinite leads. While the potential in the leads are held fixed, the 
potential in the device can spatially vary. 



A. Tight binding Hamiltonian for a one dimensional device 



Consider a system described by a set of one dimensional grid / lattice points with uniform spacing a. Further 
assume that only nearest neighbor grid points are coupled. A spatially uniform system with a constant potential has 
the Hamiltonian: 



(E-H)V = Q 



• • • 

-t E-e -t 

-t E -e -t 

-t E-e -t 

• • • 
• • • 



V • J 



(31) 



- + {E - e)* g - M q+1 = 0, 



(32) 



where, E is the energy and ^> q is the wave function at grid point q. The Hamiltonian matrix is tridiagonal because 
of nearest neighbor interaction. The diagonal and off-diagonal elements of the Hamiltonian e and t represent the 
potential and interaction between nearest neighbor grid points q and q + 1 respectively. 
The solution of eqn. (|32|) can be easily verified using Bloch theorem to be, 



E = e + 2tcos{ka) 



and the group velocity is, 



* g = e lkqa 



IDE 2at . 



(33) 
(34) 



(35) 



The uniform tight binding Hamiltonian in eqn. f) 32 j) can be extended to a general nearest neighbor tight binding 
Hamiltonian given by, 



t q , q -!^g-i +(E- e q )^ q + t M+1 ^ q+l = 0, 



(36) 



where, e q is the on-site potential at grid point q and t q q+ i is the Hamiltonian element connecting grid points q and 
q + 1. = t q q+1 because the Hamiltonian is Hermitian. 



In the special case of the discretized Schrodinger equation on a uniform grid, 

h 2 , „ h 2 



t = t 



9,9+1 tg+l,g 



(37) 



where a is the grid spacing and V q is the electrostatic potential at grid point q. 
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B. Eliminating the Left and Right semi-infinite leads 

A typical nanodevice can be conceptually divided into three regions (Fig. [6]): 

• Left semi-infinite lead (L) with a constant potential e; 

• Device (D) with and arbitrary potential, and 

• Right semi-infinite lead (R) with a constant potential e r 

The potenial of the left (right) lead ei (e r ) and the hopping parameter ti (t r ) are assumed to be constant, which 
signifies that the leads arc highly conducting and uniform. Then the Hamiltonian of the device and leads, eqn. (|36l) . 
is an infinite dimensional matrix that can be expanded as, 



• 

-tjtfj-i + (E - e,)*io - tj,d#i = (38) 

-t d ,i^io + (E - ei)*i - ii, 2 *2 = (39) 

-*1,2*1 + - e 2 )*2 - *2,3*3 = (40) 

-*n-l,n*n-l + {E- £„)*„ - t dr V rn+1 = (41) 

-t r ,d^ n + (E- e r )* r „+i - t r ^ rn+2 = (42) 

• 



where, the bullets represent the semi-infinite Left and Right leads. The subscript Im (rm) refer to grid point m in 
the Left (Right) lead. However, to find the electron density in eqn. (|26p . the wave function is only required at the 
device grid points. We will now discuss a procedure to fold the influence of the left and right semi-infinite leads into 
the device region. 

Terminating the semi-infinite Left and Right leads: The wave function in the contacts due to waves incident from the 
left lead are, 

= {e +ikin + s u e- ikin )u ln in region L (43) 
*™ = s r ie lkrn u rn in region R, (44) 

and the corresponding eigen values are [eqn. (|33l) ]. 

E-ei = 2ticos{kia) = ti{e lkia + e - tkia ) , (45) 

and similarly for the right contact, with the indices r. su and s r i are the reflection and transmission amplitudes. 
Substituting Eqs. (g21) and |g5J in eqn. $5E§ yields, 

su =t{ 1 (-ti+ti A) . (46) 

Substituting Eqs. (|4"3"|) and (|46p in eqn. (|3"9")l . we obtain, 

(E-a- t d . l e +tkia t- 1 t l4 )^! 1 + ^.2*2 = -2it d . lS in(k ia ) (47) 

Eq. (|47|) is a modification of Schrodinger's equation centered at grid point 1 of the Device (eqn. ([39]) ) to include the 
influence of the entire semi-infinite Left lead. Similarly, substituting eqn. (|44|) and E — e r = 2t r cos(k r a) in eqn. (|42p 
we get, 

Sri = t'Hd^n ■ (48) 

Now, substituting Eqs. (I44[) and (|4"B"|) in eqn. (l41~|) . we can terminate the right semi infinite region to yield, 

- tn-x^n-i + (E — e n - t d>r e ik n-% d )^ „ = . (49) 

Eq. (|4"9")) is a modification of Schrodinger's equation centered at grid point n of the Device (eqn. (JUJ)) to include the 
influence of the entire semi-infinite Right lead. 
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The influence of the semi-infinite Left and Right leads have been folded into grid points 1 and n of the device, for 
waves incident from the Left lead [Eqs. (|47|) and (|49)) ]. Now the wave function in the device due to waves incident 
from the left lead can be obtained by solving the following n dimensional matrix instead of the infinite dimensional 
matrices in eqns. (|38|) - (|42|) : 

= i L , (50) 

where, A is a square matrix of dimension n, and are n by 1 vectors. %l is the source function at (fc, E) due 
to the Left lead. Matrix A is, 

A = EI-H D - S /ead (51) 

and the non zero elements of E/ ea( 2 are, 

£w ltl = t d , l e +ikia ti\ d = T, L (52) 

S/eod„ „ = td,re +lkra t r H r> d = Efl . (53) 

£l and S/j are called the self-energies, and they represent the influence of the semi-infinite Left and Right leads on 
the Device respectively. The real part of self-energy shifts the on-site potential at grid point 1 from ei to ei +Re(Ei). 
The imaginary part of self energy multiplied by —2 is the scattering rate of electrons from grid point 1 of Device to 
Left lead (scattering rate = — 2Im[££]), in the weak coupling limit. 

In a manner identical to the derivation of eqn. (|50p . for waves incident from right contact, the wave function in the 
device (*^ } ) can be obtained by solving: 

A*™ = i R , (54) 

where is the source function due to the Right semi infinite contact. The only non zero elements of A, and ir 
are: 

A(l, 1) = E - ei - Si and A(n, n) = E — e n — S fi (55) 

A(i, i) = E-a , A{i, i + l) = -t u+1 and A(i + l,i) = -t{ i+1 (56) 

iz,(l) = —2it c i,isin(kia) and (57) 

*i?( n ) = —2itd, r sin(k r a) . (58) 



C. Electron and current densities expressed in terms of Green's functions 

The Green's function corresponding to Schrodinger's equation ([E — H]^> = 0) for the device and contacts is, 

[E-H + irj\G = I , (59) 

where r\ is an infinitesimally small positive number which pushes the poles of G to the lower half plane in complex 
energy, and H is the Hamiltonian. The Green's function of device region D with the influence of the contacts included 
is, 

AG = I, (60) 

where 

A = EI -H D -Z lead {E) (61) 



is an n dimensional matrix defined in eqns. (|55j) and (|56| . The formal derivation of eqn. (|60|) is given in the appendix 

H 

Using the definition for G in Eqs. (|50[) and (|54p , the wave function in region D due to waves incident from Left and 
Right contacts can be written as, 



^> = Gi L and (62) 

D 



*^ = Gi R . (63) 
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As il and i R are non zero only at grid points 1 and n, the full G matrix is not necessary to find the wave function in 
the device; Only the two columns G(:, 1) and G(:, n) are necessary. 

The electron density at grid point q can now be written using eqns. (|62|) and ((63]) in eqn. ([26]) as 

rig = ^ Gg^iLi^G^i^fL + G qjn i R i R G^ nA f R (64) 

k . s 

= ^2 G qA^ t d,isin 2 (kia)t^ d f L ]G\^ q + G q , n [tt d , r sin 2 (k r a)t rid f R ]G^ (65) 

k,a 

where G^ is the Hcrmitcan conjugate of the Green's function. The summation over k can be converted to an integral 
over E by, 

Using eqn. ([66]) and l^f I = 2a|t||sin(A;ja)|, eqn. (|65| becomes, 



/ S [G q AmT(E)Gi q (E) + G q , n (E)Z%(E)Gl q (E)] ■ ~ (67) 



n 9 = J 
where, 

E^(^) - 2t dii l|sm(A;;a)|t M / i (E) and (68) 

= 2t d , r ^|sm(A; r a)|t r , d /fi(£;) . (69) 

fc; and fc r at energy E 1 are determined by _E = e; — 2ticos(kio) and i? = e r — 2t r cos{k r a) (eqn. (I33p ) respectively. It 
can be seen from Eqs. ([52]), ([53]), jM]) and ((Ml): 

Ef(£) = -2Im[V L {E)]f L {E) (70) 
E»(£) = -2/m[S R (S)]/ i? (£;) . (71) 
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The electron density [eqn. I|67p] can then be written as, 



2tt 



G{E)T?£ ad {E)G^{E)\ q , q .-, (72) 



where, the non zero elements of £™ ad are 



E- dia (i?) = and (73) 

T%ad n JE) = ^r(E). (74) 

and £^ defined above in Eqs. ([7D|) and (fTTj) arc called the in-scattering self-energies due to contacts. These 
self-energies physically represent in-scattering of electrons from the semi-infinite leads to the device and so play an 
important role in determining the charge occupancy in the device. They depend on the Fcrmi-Dirac factor / occupancy 
in the contacts /z, and f R , and the strength of coupling between contacts and device, Im[T,L(E)] and Im[Ti R (E)]. 
It is easy to see using eqns. (fT2"j) and (|?4"|) that the electron density can also be written as, 

n q = f dE[LDOS L (q,E)f L {E) + LDOS R (q,E)f R (E)}, (75) 

where, LDOSL(q, E) (LDOS R (q, E)) are the density of states due to waves incident from the left (right) contact, at 
grid point q, and 

G n 1 r T J G\ n 1 

LDOS L {q,E) = qA hi .- (76) 
7r a 

LDOS R (q,E) = G ^ TrG ^ . 1 , (77) 



where, 



T L (E) = -2Im[E L (E)} (78) 
T R (E) = -2Im[X R (E)} . (79) 



Note that eqn. (|75|) is identical to eqn. (j2~l]) for scmiclassical ballistic transport. 



The current density beween grid points q and q + 1 per unit energy can be written as [eqn. (|27[) ] , 

Now, following the derivation for electron density above [eqn. I67| . it is straight forward to derive that the current 
density 

-G q+hl (E)ZT(E)G\(E) G q+hn (E)X%(E)Gl q (E)]. (81) 



The current density is given by 



Electron correlation function: More generally, we define the electron correlation function G n , which is the solution 
to 

AG n = E}» d & . (83) 
Noting that G = A' 1 , it is easy to obtain eqn. (J72J) G n = GS™ Qd G t . 
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The expressions for electron [eqn. (|72[) ] and current [eqn. (|52"j) ] densities at energy E, in the phase coherent case at 
finite applied biases can now be written as, 

n q (E) = 2-^-1 , (84) 

"W 1 ^ = 2^ 2 ^ [G ^ +l{E) - G ^ m - (85) 

That is, the diagonal and first off-diagonal elements of G n are related to the electron and current densities respectively. 
Note that these equations are entirely equivalent to eqns. (|26|) and (|27|) appearing in the Landauer-Buttiker approach. 



Hole correlation function: In the absence of phase breaking scattering, the Green's function (G) and the electron 
correlation function (G n ) are sufficient for device modeling. Scattering introduces the need for the hole correlation 
function (G p ), whose role will become clearer in section IV Dl While the less-than Green's function is directly pro- 
portional to the density of occupied states, the hole correlation function is proportional to the density of unoccupied 
states. 

The density of unoccupied states at grid point q is also obtained by applying the Landauer-Buttiker formalism. For 
this, we simply replace the probability of finding an occupied state in the contact /l,r by the probability of finding 
an unoccupied state in the contact, 1 — /l,_r, in eqn. (|2"6"|) . Then following the derivation leading to eqn. (|2"6"|) . we 
obtain: 



dE r 
2^ 
dE 



1 



h i = 2 / T^G(E)^ d G\E)\ 



2tt 



(86) 
(87) 



where, the only non zero elements of are, 



(E) 



:° L ut (E) = 2t dtl -sin(ha)ti, d (l ~ f L ) = 
= 2t d , r \sin{k r a)t rA {\ - f R ) 



-2Im[X L (E)][l-f L (E)} 
-2Im\Z R (E)][l- f R (E)]. 



(88) 
(89) 



Akin to Eqs. (|83|) and (|84|). the density of unoccupied states at energy E, at 
diagonal elements of G p , 


grid point q can be expressed as the 


hq ( E )=2 2 ^E), 


(90) 


where, G p is in general given by, 




AG p = S^adG 1 • 


(91) 



D. Electron-phonon scattering 

In section IV C\ we defined the retarded, in-scattering and out-scattering self-energies in the device arising from 
coupling of the device to the external contacts. The self-energy S™ represents in-scattering of electrons (in-scattering 
rate) from the semi-infinite left contact to the device, assuming that grid point 1 of the device is empty. A similar 
statement applies to • The in-scattering self-energy of the contacts depend on their Fermi distribution functions 
and surface density of states. 

A second source for in-scattering to grid point q and energy E [(q, E)] is electron-phonon interaction. The self- 
energy at (q, E) has two terms corresponding to in-scattering from (q, E + ftio p } lonon ) and (q, E — hio p honon), as shown 
in Fig. [71 Intuitively, the in-scattering self-energy (in-scattering rate) at (a;, E) should depend on the Bose factor 
for phonon occupancy, the deformation potential for electron-phonon scattering and the availability of electrons at 
energies E + ftwphonon and E — hiUphonon- It follows rigorously, within the Born approximation that the in-scattering 



16 



less-than / in- scattering self energies Z D 
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Left Lead Device Right Lead 



FIG. 7: Pictorial representation of the meaning the two in-scattering self-energies that appear in this chapter. S^(£J) and 
Y^(E) are the self-energies of the leads, and are non-zero only at the first and last Device grid points. Epf lonon ln (i?) is the 
self-energy due to electron-phonon interaction, and is non-zero at all Device grid points. 



self energy into a fully empty state at energy E and grid point q is }Mah87l | 

oiioi I, on on on on 01 ion )] ■ (92) 

v 

represents the electron-phonon scattering strength at grid point q. The first term of eqn. f52")) represents in- 
scattcring to E from E — hw v honon (phonon absorbtion) . n B is the Bose distribution function for phonons of energy 
fujJphonon and G™(E — hjjphonon) is the electron density at E — hw v honon- The first and second terms of eqn. (|9"2"|) 
represents in-scattering of electrons from E — hu> p honon (phonon absorbtion) and E + hujphonon (phonon emission) to 
E respectively. The in-scattering rate at grid point q is given by 

In-scattering rate at grid point q: — — - = T, q n (E) , (93) 



where £™ is the sum of all in-scattering self energies at grid point q. 

The out-scattering self-energy E™ 4 in eqn. (|88|) represents out-scattering of electrons from grid point 1 in the device 
to the semi-infinite left contact, assuming that grid point 1 of the device was occupied. The out-scattering self-energy 
due to the left contact depends on the probability of finding an unoccupied state in the left contact 1 — /l and the 
surface density of states of the left contact. A similar statement applies to E™*. A second source for out-scattering of 
electrons from an occupied state at (g, E) is electron-phonon interaction, which leads to scattering to (g, E + hujphonon) 
and (g, E — Hu> phonon) as represented in Fig. [S] Intuitively, the out-scattering self-energy (out-scattering rate) at (q, E) 
should depend on the Bose factor for phonon occupancy, the deformation potential for electron-phonon scattering 
and the availability of unoccupied states at enegies E + ftw p honon and E — hujphonon ■ It foll ows rigo rously, within the 
Born approximation that the out-scattering self energy from a fully filled state at (q, E) is |Mah87l | 

Z°p U honon q ( E ) = D Q [i 71 B phonon) + l)G P q (E - hu phonon ) + n B (hu phonon )G P q (E + frw phonon )] . (94) 

In the above equation, G P (E ~ hai p } lonon ) and G P (E + hui p i lonon ) are the densities of unoccupied states at E—hu> p honon 
and E + fuo v h onon . So the first and second terms of eqn. (|94[) represents out-scattering of electrons from E to 
E + hujphonon (phonon emission) and E — hujphonon (phonon absorbtion) respectively. The out-scattering rate at grid 
point q is given by 

Out-scattering rate at grid point q: — = E°"*(i?) , (95) 

where E° u * is the sum of all out-scattering self energies at grid point q. 

We now discuss how the in-scattering self-energy due to electron-phonon scattering affects the expression for electron 
density. The electron density at grid point q in the phase-coherent case (eqn. (|72[l) is the sum of two terms, 

« f dE 
n n = 2 — 



J £ [Gq,i(mT(E)Gl q (E) + G qA (E)T^(E)Gl q {E)] ■ ~ . (96) 
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FIG. 8: Pictorial representation of the meaning the two out-scattering self-energies that appear in this chapter. T,i ea d° ut (E) 
is self-energy due to leads, which is non-zero only at the first and last Device grid points 'Ephonon° ut {E) is self-energy due to 
electron-phonon interaction, which is non-zero at all Device grid points. 



The first term represents in-scattering of electrons from the left contact [(E™(i5)), which is propagated to grid point 
q via the term G q ^i(E)G[ q (E). The interpretation of the second term is similar except that it involves the right 
contact. In the presence of electron-phonon interaction, the in-scattering functions Sp /lonon is non zero at all grid 
points. As a result, an electron can scatter from (q',E') to {q',E) and then propagate to grid point {q,E) via the 
term G q ^ q '{E)G\ g (E). The expression for the electron density can be generalized to include such terms, 



2/^ 

2tt 



g [G(E)^ ad (E)GHE) + G{E)^ honon {E)G\E)] \ M ■ i 



(97) 



dE 



1 



where the third term corresponds to propagation of electrons from grid point q' to q after a scattering event at q' as 
shown in Fig. ([9]). The in-scattering self energies due to phonon scattering are given by eqn. (|92[) . More generally, 
G" is given by, 



G n (E) 

[E - H - Y,(E)}G n (E) 



G(E)T, m {E)G t (E) 
Y, in (E)G^ (E), 



(99) 



where S ln is the sum of self-energies due to leads, electron-phonon interaction and all other processes. The reader 
can compare the above two equations to Eqs. (|73[) and (|83p . which are valid in the phase coherent limit. 
The density of unoccupied states can be written in a manner identical to eqn. (|97[) as, 



2/^ 

2tt 



G q , 1 (E)^° L ut (E)Gl q (E) +G q , 1 (E)^ t (E)G{ q (E) +Y,G q ^{E)^ Phonon {E)Gl, ^ q (E) 
= 2 y H [G{E)^ d {E)G\E) +G(E)Z° P t onon (E)GHE)] \ M ■ ~ 

More generally, the G p matrix is given by, 



(100) 



G P (E) = G{E)Y, out {E)G Js {E) 
[E - H -Z(E)]G P (E) = Y, out (E)G*(E), 



(101) 
(102) 



where S otl * is the sum of self-energies due to leads, electron-phonon interaction and all other processes. 
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Three contributions to electron density in eqn. 97 
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FIG. 9: Contributions to electron density from leads and electron-phonon scattering. 
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FIG. 10: (Top) Examples of the layered structures: carbon nanotube, DNA molecule, MOSFET. (Bottom) Scheme of simulation 
domains in the layered structure: device, left and right leads. 



Note that in general G n and G p are full matrices, the diagonal elements of which correspond to density of occupied 
and unoccupied states respectively, and the first off diagonal elements of G n and G p correspond to the current density. 
The Green's function G in the device region is obtained by solving, 

[E — H — Y<i ea d{E) — T,phonon{E)\ G = I , (103) 

which is similar to eqn. (|60p for the phase coherent case, except for the additional self-energy due to phonon scattering. 



VI. NON-EQUILIBRIUM GREEN'S FUNCTION EQUATIONS FOR LAYERED STRUCTURES 

The previous section dealt with a simple one dimensional Hamiltonian. In this section, we will present the NEGF 
equations for a family of more realistic structures called layered structures. A layer can be considered to be a 
generalization of a single grid points / orbital (Fig. to a set of grid point / orbitals (Fig. [TO]) . Three examples of 
layered structures are shown in Fig. [TO] The left most figure is a carbon nanotube, the middle figure is a DNA strand 
and the right figure is a MOSFET. A layer consists of the atoms / grid points between the dashed lines in Fig. [TO] 

A common approximation used to describe the Hamiltonian of layered structures consists of interaction only between 
nearest neighbor layers. That is, each layer q interacts only with itself and its nearest neighbor layers q — 1 and q + 1. 
Then, the single particle Hamiltonian of the layered structure is a block tridiagonal matrix, where diagonal blocks H q 
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represent the Hamiltonian of layer q and off-diagonal blocks T q ^ q+ i represent interaction between layers q and q + 1: 

\ 



• Hi T 12 

t\ 2 h 2 t 2 . 3 



H 



T 



n-2,n-l 



T 



t 



V 



n — l.n 



T n -i. 



(104) 



The Green's function equation in the presence of electron-phonon scattering is 

[EI — H — ^Phonon]G = I , 



(105) 



is the self-energy due to electron-phonon scattering. We partition the layered structure into Left lead, 
Device and Right lead as shown in Fig. [TUJ The Device corresponds to the region where we solve for the nonequilibrium 
electron density, and the leads are the highly conducting regions connected to the nanodcvice. While the Device region, 
where we seek the non equilibrium density consists of only n layers, the matrix equation corresponding to cqn. (|105[) 
is infinite dimensional due to the semi-infinite leads. We will now show how the influence of the semi-infinite leads 
can be folded into the Device region. In a manner akin to the previous section, the influence of the semi-infinite 
leads is to affect layers 1 and n of the Device region. An important difference is that the derivation here includes 
electron-phonon scattering and does not assume a flat potential in the semi-infinite leads. We first define 



A' = [EI - H - £ 



Phononl 



Then, eqn. (|105|) can be written as, 



A'G = I 
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(106) 



(107) 



where, 



A i 



LL 



4' 



V 

( A' rl 
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4' 
1 n 



A' 



corresponds to the left semi-infinite lead, 



(108) 



RR 



— T r i 
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corresponds to the right semi-infinite lead, and 



(109) 
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corresponds to the Device region. (110) 
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(111) 
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are the coupling between the Left and Right leads and Device respectively. Note that A' DL = A' LD ^ , A' DR = A' RD \ 
and A' LD and A' DL {A' RD , and A' DR ) are sparse matrices. Their only non-zero entry represents coupling between the 
Left (Right) lead and Device. O represents zero matrices. From eqn. (|107p . we have: 

G LD = -A^A' LD G DD (112) 

Grd = — A' rr A' rd Gdd (113) 

A'dlGld + A' dd Gdd + A dr Grd = I ■ (114) 

Substituting Eqs. (|112p and Q113P in eqn. (|114[) . we obtain a matrix equation with dimension corresponding to total 
number of grid points / orbitals in the n device layers , 

[A'dd ~ A' DL A' LL A' ld — A' dr A' rr A' rd ]Gdd = I ■ (115) 

The second and third terms of eqn. (|115|) are self energies due to coupling of the Device region to Left and Right 
leads respectively. 

The Green's functions of the isolated semi-infinite leads by definition are, 

A' LL g L = I and A' RR g R = I . (116) 

The surface Green's function of the Left and Right leads are the Green's function elements corresponding to the 
edge layers — 1 and n + 1 respectively, 

ffi-L-i = A 'll 1a ^d 9R n+1 _ n+1 = A RRi i . (117) 
Eq. (|115[) can now be rewritten in a form very similar to eqn. (|60[) . 

[EI — H — Yiphonon — ^>lead]GDD = I (H8) 

where, 

Sjeodi,i = Tdl9l_ 1i _ 1 Tld = (119) 

^leadn,n = T DR g Rn + 1 n + 1 T RD = Y, R . (120) 

All other elements of S; ea£ j are zero. and are self-energies due to the Left and Right leads respectively, and 
Tdl =T£ d and T DR = T RD . Finally, defining, 

Add = EI — H — Y^phonon — ^lead , (121) 



eqn. (|118[) can be written as 

A DD G DD =I. (122) 



The main information needed to solve eqn. (|118p are the surface Green's functions of gL and g R . We will disucss 
two methods to obtain this surface Green's functions for a constant potential in the Left and Right leads. When the 
potential does not vary, A' LL and A' RR are semi-infinite periodic matrices with all diagonal / off-diagonal blocks being 
equal: 

A n =A l2 = A l3 = ... = Ai (123) 
T n =T l2 = Ti 3 = ...=Ti . (124) 

gh_ 1 1 is obtained by solving the matrix quadratic equation: 

[Ai-TUl_^Ti]9l_^=I ■ (125) 
This equation can be solved iteratively by, 

[JU - T}g L <Z?>Tt]g L <Z> = I , (126) 

where, the superscript of gL represents the iteration number. Note that the solution to eqn. (|125|) is analytic when 
the dimension of Ai is one. A second simpler solution to obtain g^ involves transforming to an eigen mode basis using 
an unitary transformation (S), such that 

S^AiS = A M and S^TiS = T, d , (127) 
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where, both Aid and are diagonal matrices. The surface Green's function in this new basis is simply a diagonal 
matrix, whose elements are obtained by solving the scalar quadratic version of cqn. (|125[) . The Green's function in 
the original basis (in which Ai is not diagonal) can be obtained using the inverse unitary transformation. 

Electron (G n ) and Hole (G p ) Green's Function: 

The electron density is equal to [see the discussion of electron density in section IV C| , 

G n (r,f,E) , x 

n(f,E) = 2 y - > . (128) 

The governing equation for G n is 

[EI - H - T,phonon]G n = T, l p honon G f , (129) 

where is the hermitian conjugate Green's function and EpV on is the in-scattering self-energy due to phonon 
scattering. The dimension of Eq. (| 1 29[) is essentially infinite due to the semi-infinite Left and Right leads. It can 
however be converted to a finite dimensional matrix with dimension equal to the number of grid points / orbitals 
corresponding to the n device layers. In a manner identical to the derivation of eqn. Q118P for the retarded Green's 
function, it can be shown that the role of the leads can be folded into layers 1 and n to yield, 

[EI — H — ^Phonon — ^lead]G DD = ^DD^DD ' (130) 



A D DG n DD - ??£ D G DD , (131) 

where Ado has been defined in eqn. (|12ip . The self-energy S^^, has contributions due to both electron-phonon 
interaction and leads, 

(132) 
(133) 
(134) 

The in-scattering self-energies due to the leads, and have forms very similar to eqns. (|70|) and (|7Tj) of 
section IV CI 

Ef(E) = -21m[E L (E)]f L (E)=r L (E)f L (E) (135) 
= -2 Im[V R (E)]f R (E) = T R (E)f R (E) , (136) 

where, 

T L (E) = -2Im[S i ( J B)] (137) 
T R (E) = -2 Im[E R (E)] . (138) 

/l and f R are the distribution functions in the Left and Right leads respectively (Fermi factors at equilibrium). The 
self-energies T,l(E) and T, R (E) have been defined in eqn. (|119[) and (|120p . 

The diagonal elements of the hole correlation functions G P (E) represents the density of unoccupied states, 

Mr -, £) = 2 5!MS. < 139) 

The governing equation for G P (E) is, 

[EI — H — ^phonon]G p = Y, P u l * onon G' ! (140) 

or 

A DD G P DD = Y% D G DD (141) 

•sr\out . \-\out I •spout f'\A r )\ 

^DD^ ~ ^Phonon^ \ 14Z J 

spout spout I snout /i ao\ 

^DD^ ~ ^Phonon nn ~r V 1 ^) 

ES-* M = S™ 4 „ „„ , where, i=2, 3, 4, ... n-1 (144) 



spin 

lj DD 1A 


spin 

Phonon i i 






spin 


spin 

Phonon nn 






spin 

DD <i, 9 


spin 

Phonon qq 


, where, q=2, 3, 4, .. 


. n-1 
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where ^°pjl onon is the out-scattering self-energy due to phonon scattering. The out-scattering self-energies due to 
the leads, E£ and Eg*' have forms very similar to eqns. ([55)) and (1591) 

Ef = -2 Im[E L (£)](l - f L (E)) = T L (l ~ h{E)) (145) 
= -2 Im[EH(S)](l - /«(£)) = r„(l - /«(£)) , (146) 

where 1 — Jl(E) (1 — Jr(E)) is the probability of finding an unoccupied state in the left (right) contact at energy E. 
Eqns. (fT3"Tj) and (fiTTj) for G n and G p can be written as. 

GdD = ^DD^DD^DD = GdD^dD^DD (147) 

G P dd = ^dd^Sd^Ld = Gdd^dd&dd ■ (148) 

While these equations appear often in literature, we do not suggest using them to compute the diagonal elements of 
G n and G p of layered structures. This is because their use requires knowledge of the entire Gdd matrix when 
is non zero at all grid points. Computation of the entire Gdd amounts to inversion of Add- Matrix inversion is 
computationally expensive and scales as N 2 - 7 where N is the dimension of Add- The diagonal elements of G n and G p 
of layered structures can be computed more efficiently without calculating the entire Gdd matrix using the algorithm 
discussed in appendix [Cl 

Current Density: We will now present some expression for current density commonly used in literature. The current 
flowing between layers q and q + 1 is (eqn. (|55"|) of section IV CP : 



ie f dE 

Jq->q+l = T-2 / -^r Tr [Tq,q+lGq +ltq {E) - T q+1 , q G qq+1 (E)] 



(149) 



Eq. (|149[) frequently appears in the literature in other useful forms that are derived below. Expanding both terms of 
eqn. (I149|) using eqn. (|C10[) of appendix [Bl we get, 

Jl = -2 J gTr([T iD G 14 (£;)T DL5 E o j£;)+T LZ5 Gr a (£;)T DL . g [ o o (£;)] 

-[TM^^J^T^G^^ + TMS^^r^GJ,!^)]) (150) 
= f 2 / fir ([G,,^) - G[ ^BjlT^s^^mD - q^EjT^b^^B) - ffi 00 (B)]T M ) (151) 
Using the relationships, 

Si" - Tdl91 Tld (152) 
-iT L = T DL [g Lo o -g{ o jT LD , (153) 

Eq. (|151|) can be written as, 

Jl = ^J ^Tr(i[G ltl (E) - G\ A (E)]Xt(E) G^)^)) (154) 

Eqns. (|149p and (|154[) are bo th general expression for current density valid in the presence of electron-phonon 
scattering in the device |Mei90l | . The advantage of using eqn. (|149|) is that the current density can be calculated 
at every layer in the device. This expression is useful in understanding how the current density is energetically 
redistributed along the length of the device as a result of scattering. 

In the phase coherent limit, T,phonon = and the only non zero self-energies are in layers 1 and n due to the 
contacts. We define matrices, Y ^ and Tr, which consist of n blocks corresponding to the n Device layers (dimension 
of Add matrix) and with the following non zero elements: 

f L | 11= r L , f R \ nn =T R . (155) 

Now left multiplying eqn. (|122|) by G^ and right multiplying the Hermitian conjugate of eqn. (|122p by G, and 
subtracting the resulting two equations, we have, 

G-G f = G t (E-£ t )G , (156) 
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where S is the total self-energy due to phonons and the leads. The Hcrmitcan conjugate Green's functions and 
self-energies are and £T . In the absence of phonon scattering, the self energies only have components due to the 
leads and so cqn. ()156|) can be written as. 

i[G - G f ] = G f (f L + f R )G , where (157) 

eqns. (|T37)) and (fl"38)) have been used. It also follows from Eqs. ([135]) . (fl"36)) and (fl"47)) that 

G n = G (f L f L + T R f R )G^ . (158) 

Now using eqns. (|157[) and (|158[) in eqn. ()154[> . the current in the phase coherent limit is, 

Jl = ^J ^T(E)[f L (E) f R (E)} . (159) 

The total transmission at energy E is identified from eqn. (| 1 59[) to be 

T(E) = Tr[f L (E)G(E)f R (E)G^(E)} . (160) 

Note that to compute the total transmission using eqn. (|160| . only the elements of G connecting layers 1 and n are 
required because Tl and T R are non zero only in layers 1 and n respectively. 



A. Crib Sheet 



The algorithmic flow in modeling nanodevices using the non equilibrium Green's function consists of the following 
steps (Fig. [TTj) . We first find a guess for the electrostatic potential V(r) and calculate the self energies due to the 
contacts (eqns. (|175|) - f| 185[) ^) . The self-energies due to electron-phonon scattering are set to zero. The non equilibrium 
Green's function equations for G, G n and G p (eqs. (|171[) - cqs (|174jl ) are then solved. Following this, the self energies 
due to electron phonon scattering and contacts (eqns. (|175[) - (|185[) )are calculated. As the equations governing 
the Green's functions depend on the self energies, we iteratively solve for the Green's function and self-energies, as 
indicated by the inner loop of Fig. QT] Then, the electron density (diagonal elements of G n ) is used in Poisson's 
equation to obtain a new potential profile. We use this updated electrostatic potential profile as an input to solve for 
updated non equilibrium Green's functions, and continue the above process iteratively until convergence is achieved 
(outer loop of Fig. [Tl]) . A number of equations that are repeatedly used in nanodevice modeling are listed below. 

Physical Quantities: 



Scattering Rate: 



-(E) 



= -2 Im[Z(E)] = T(E) 



(161) 
(162) 



Density of States at (r, E): N(r, E) = Im(G(r, r, E)) 



(163) 



Use recursive algorithm to calculate DOS (Do not invert A). 



/dE 
—G n (r,r,E) (164) 

Use recursive algorithm to calculate n (Do not use G n = GE m G^). 



Unoccupied Density at 



h{r)=2 j^GP{r,r,E) (165) 
Use recursive algorithm to calculate h (Do not use G n = GY, in G^). 
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Guess for potential profile 



1 



Solve for contact self-energies 
Solve Non-Equilibrium Green's Function Equations (G, G n , G p ) 
Input: Device Geometry, Potential profile 
Output: Electron density, Current, G, G n , G p 



1 



Equations governing scattering 
Input: Electron-phonon deformation potentials 
Output: Scattering rate, Self-energies E, E ln , E 0Ljt 



Poisson's equation 
Input: Charge & doping densities, boundary conditions 
Output: Potential profile 



1 



No 



Converged? 



1 



Yes 



Current & electron densities, Scattering rates 



iterate 



Number 
crunching step 

Solve at 
each energy 



FIG. 11: Flowchart of a typical simulation involved in modeling of a nanodevice. 



Current density flowing between layers q and q + 1 (valid with scattering in device): 

ie f rlE 

J q ^ q+1 = j2j — Tr [T qtq+1 G n q+1>g (E) - T g+M G^ +1 (£)] (166) 
Current density flowing from the Left lead into layer 1 of Device (valid with scattering in device): 

•h = \lj ^Tr{i[G hl (E)-Gl 1 (E)]Ef(E)-Gl 1 (E)T L (E)} (167) 

= I 2 ! ^{Glii^TW-Gl^E^iE)} (168) 
Current density flowing from the Left lead into layer 1 of Device (valid only in the phase coherent limit): 

Jl = J? I ^ T(E) [h{E) f * m > (169) 
where the total transmission from the Left to Right lead at energy E is given by 

T(E) = Tr[Y L {E)G(E)t R (E)G\E)]. (170) 
Only elements of G connecting layers 1 and n are necessary. 

Equations Solved: 

Green's Function: [EI — H — H]G{E) = 7 — » AG = I (171) 

Hermitean conjugate Green's Function: G f (E)[EI — (172) 

electron correlation Function: [EI - H - Z}G n {E) = T, in {E)G^{E) -> AG n = T, in G 1 (173) 
hole correlation Function: [EI - H - T,]G P {E) = S ollf \E)G\E) -> AG P = Y, out G^ (174) 
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X a (E) = ^ead(E) + ^P honon (E), where a G in, out (175) 

Ef*^ - ?, a L {E)=T DL gl o T LD (176) 

Ew„,„ = SS(^) = T D R9% n+ln+ T RD (177) 

Sw M =0 V i?l,n (178) 

r(S) = -2Im[£(£T)] (179) 

T L (E) = -2Im{E L {E)} (180) 

T R (E) = -2Im[E R (E)} (181) 

E™(£) = T L (E)f L (E) (182) 

= T R (E)f R (E) (183) 

sr 4 (S) = r L [l-/ L (£)] (184) 

V R \E) = T R [l-f R (E)} (185) 



The diagonal and nearest neighbor off-diagonal elements of G and G n are computed repeatedly as they correspond 
to physical quantities such as the density of states, electron density and current. Non local scattering mechanisms, 
which requires calculation of further off-diagonal elements, are not discussed here. 

Useful Relationships: 



z[G-G f ] = G p + G n = DOS (186) 

i[£-£ f ] = S ou * + S ira = r (187) 

G-G< = G[£ f -E]G 1 (188) 

= Gtpt-EJG (189) 

= -iGrG f = -?:G f rG (190) 

G nt = G™ (191) 

G pt = G p (192) 



VII. APPLICATION TO A BALLISTIC NANOTRANSISTOR 



Quantum mechanics is playing an increasingly important role in modeling transistors with channel lengths in the ten 
nanom eter regime for many reasons: (i) Tunneling from gate to channel and source to drain determine the off current 
Svi02| . (ii) Ballistic flow of electrons in the channel is important as the channel length becomes comparable to the 



electron mean free path, (iii) Classically, the electron distribution in the inversion layer is a sheet charge at the Si-SiC>2 
interface. Quantum mechanically, the inversion layer charge is distributed over a few nanometers perpendicular to 
Si-SiG"2 interface due to quantum confinement. Methods based on the drift-diffusion and Boltzmann equations do 
not apriori capture the above mentioned quantum mechanical features. 

In this section, we will disc uss the equations involved in the two dimensional modeling of nanotransistors within 
the effective mass frame work Ren03| . The quantum mechanical and semiclassical results are compared to illustrate 



their differences. The discretized matrix equations that we solve are discussed in section IVII A[ The application of 
the quantum mechanical method to illustrate the role of polysilicon gate depiction, and the slopes of Id versus gate 
(V g ) and drain (Vd) voltages are discussed in section IVIIBI 



Table 1. List of Abbreviations: Length Scales 




oxide thickness 


hp 


polysilicon gate thickness in x-direction 




boundary of substrate region in x-direction 


Ly 


Poisson's and NEGF equations are solved from 
-Ly/2 to +L y /2 


L g 


length of polysilicon gate region in y-direction 
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A. Discretized equations 




• • • 



-10 1 



semi-infinite 
boundary 



q q+1 



• • • 



Z Z 



+L 



B 



semi-infinite 
boundary 



FIG. 12: The equations are solved in a two dimensional non uniform spatial grid, with semi-infinite boundaries as shown. Each 
column q corresponds to the diagonal blocks of the Green's function equations. 

We consider N b independent valleys for electrons within the effective mass approximation. The Hamiltonian of 
valley b is 

rr ,^ h 2 r d / 1 d\ d ( 1 d\ d ( 1 d Y » « 

Hb ^ = -T [ T X {^ x T X ) + Ty{ri>d-y) + d- Z ) ] + V ^ (193) 

where (m^,m^,m^) are components of effective mass in valley b. The equation of motion for the Green's function 
(G) and electron correlation function (G n ) are 

[E - H bl (ri )] G bl , 6a (fi ,r 2 ,E)- J dr E bl , b , (n , r , E)G b ,,b 2 (r , f 2 ,E) = S bl , b J(ri - f 2 ) (194) 

and 

Mp) (st. ^ t?\ _ / ^v. ..(st. s? r?\r> n (p) 



[E-H^G^^in^E)- I d?X bub ,(? 1 ,?,E)G^ 2 {r,f 2 ,E) = 

dr (n , r, E)Gl M (r, r 2 ,E). (195) 



The coordinate spans only the device in Eqs. (|194p and (|195p . The influence of the semi-infinite source (S), drain (D) 
and polysilicon gate (P) leads, and electron-phonon interaction are included via self-energy terms f and j]™(° M *) as 
discussed in section PVTl The contact self-energies are diagonal in the band index, E" &9 c = E" c S bl , b2 (C represents 
contacts) . 

The electrostatic potential varies in the (x, y) plane of Fig. 1121 and the system is translationally invariant along 
the z-axis. So, all quantities A(ri,r 2l E) depend only on the difference coordinate z\ — z 2 . Using the relation 

A(n,?2,E) = J ^e lk ^-^A( Xl , yi ,x 2 ,y 2 ,k Zl E) , (196) 
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the equations of motion for G and G n ^ simplify to 



and 



[E - ^ - H b {h)]G b {ri,r 2 ,k z ,E) - / dFH b (n,f, k z ,E)G b (r,f 2 , k z ,E) = 8(n - ffe) (197) 



[£-^-ff6(ri)]G 6 (fi,f 2 ,fcz,£)- / df^ b (r 1 ,r,k Zl E)G n b (v \r,r 2 ,k z ,E) 

ZJTl-y 



drZ^^irur, k z ,E)G\{r, r 2 ,k Z7 E), (198) 

where Z b = Z bjb , and r — > (a;, y) for the remainder of this section. 
Eqs. (|197[) and (|198[) can be written in matrix form as, 

A'G = A and (199) 
and A'G n(p) = s» l (°«*)G*t . (200) 

The self-energies due to S, D and P are non zero only along the lines y = ys = yi, y = Ud = J/jv„ and £ = —(Lp + t ox ) 
respectively of Fig. [TSl The A' matrix is ordered such that all grid points at a y-coordinate (layer) correspond to a 
diagonal block of dimension N x , and there are N y such blocks. In the notation adopted A' 3i - 2 (i, i') is the off-diagonal 
entry corresponding to grid points (x^y^) and (x^ , yj 2 ) . The non zero elements of the diagonal blocks of the A' 
matrix are 

A' jd (i, i) = E' - V itj - T jd (i + - T jd (i - l,i) - T j+ltj (i, i) - Tj_ ltj (i, i) 

-'E s (x i ,x i )S jtl - Ei?(a; < ,x i )5j- ) jv s - T, P {y j ,y j )5 itl - £(x l; x 4 , ?/,) (201) 
± = T jtj (i± - T lS {x i±1 ,x i )5 jA - Z D (x l±1 ,Xi)5 jtNy 

-^{x i ±i,y 3 ;x i ,y J ) (202) 
(*'*') = ~ ^(xj,^')^-,^, for i' ^ i, i± 1 , (203) 

where E 1 = E — h 2 kl/2m z and Vij = V(xi, yj). The off-diagonal blocks are 

^■±i,j(M) = 2 i±i,j(i,*) _s p(fj)2/i±i) <J i,i 
A^CM') = O.for/^j, i±l, (204) 

and the non zero elements of the T matrix are 

2jj(i±l,i) = , , (205) 

2m^ x x i+ i - Xi-i \Xi±i ~ Xi\ 

T j± ij(i,i) = A: , r, (206) 

2m ± » y j+1 - yj-i \y j±1 - y 3 \ 

where m ±x = m,±1,J 2 +m ' ,J and m ±v = m, ' ;±1 2 +miJ , Non zero elements of Sp(j/j, y^), where j' ^ j,j ± 1 are neglected 
to ensure that A' is block tridiagonal; The algorithm to calculate G and G n relies on the block tridiagonal form of A 1 . 
The A appearing in eqn. (|199[) corresponds to the delta function in eqn. (|197[) . A is a diagonal matrix whose elements 
are given by 

= 7 1? \ ■ ( 207 ) 

(x i+ i - Xi- 1 ){y i+1 - yi-t) 



B. Results 



We will now discuss quantum mechanical aspects of transport in a two dimensional ballistic nanotransistor. We 
do so by comparing the current-voltage characteristics from our quantum and drift-diffusion simulations as shown in 
Fig. [T3l The important features are a higher off-current, threshold voltage shift, smaller subthreshold slope and much 
higher on-current, in the quantum case. 
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FIG. 13: Plot ol drain current versus gate voltage from the quantum mechanical calculations and Medici, at Vd = IV. At small 
gate voltages, the drain current from Medici are comparable to the 'Ql fiat band' results. The drain current from 'Ql q-poly' 
is however significantly different at large gate voltages. 



The change in threshold voltage results from the very different quantum mechanically calculated potential profile 
in the polysilicon gate. The quantum mechanical band bending is opposite to the drift-diffusion case as shown in 
Fig. Q3] In the quantum case, the conduction band at the polysilicon-Si02 is lower by approximately 130 meV. The 
physical reason for the qualitatively different quantum potential profile arises due to the tiny quantum mechanical 
probability density for electron occupancy close to the barrier. As a consequence, the electronic charge density is 
smaller than the uniform background doping density, near the SiC>2 barrier in the polysilicon region. This causes the 
conduction band in the polysilicon gate to bend in a direction opposite to that computed classically. 




FIG. 14: Potential profile at the y=0 slice of MIT25, calculated using quantum and drift-diffusion methods by assuming flat 
band in the polysilicon gate. 

The quantum mechanically calculated Id versus V g with and with out the quantum mechanical band bending in the 
gate region is shown in Fig. 1151 The gate voltage shift is approximately equal to the band bending in the polysilicon 
gate. A shift in Id(V g ) from the flat band result by the equilibrium ID built-in potential does a reasonable job 
(triangles in Fig. I15p of reproducing the results obtained by quantum mechanical treatment of the polysilicon region. 
This is especially true at small gate voltages and becomes progressively poorer with increase in gate voltage. 

The subthreshold slope d[log(Id)]/dV g is smaller in the quantum case compared to the drift-diffusion case. To 
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FIG. 15: Drain current versus gate voltage for Vd = 1 V. Quantum mechanical treatment of the polysilicon gate results in much 
higher current (solid line). The triangles correspond to the Id(V g ) calculated using a flat band in the polysilicon region shifted 
by the equilibrium built-in potential of 130 meV in the polysilicon region. 



understand the reason for this, we first note that the subthreshold current resulting from the simple intuitive expression 

I = I q0 (208) 

matches the quantum result quite accurately. I q o is a prefactor chosen to match the current at V g = 0. E r \ is 
the energy of the source injection barrier due to the lowest resonant level in the channel, which varies with gate 
bias. The higher resonant levels do not carry appreciable current. The difference in slope between the classical and 
quantum results for Id{V g ) is because of the slower variation of E r \ in comparison to the drift-diffusion barrier height 
(Ei,(classical)) as a function of V g (Fig. [TB)) . We also find that the decrease of E r \ with increase in gate voltage is 
slower than the barrier height determined from the quantum potential profiles. The slower variation of E r i arises due 
quantum confinement in the triangular well in the channel that becomes progressively narrower with increase in gate 
voltage (Fig. [T6|) . This change in confinement is a non issue in the classical case. 




FIG. 16: Left: Location of the first resonant level E r i versus gate voltage and the classical source injection barrier Eb(classical). 
Note that E r i decreases slower than Eb(classical) with gate voltage due to narrowing of channel potential well. Right: Narrowing 
of the triangular well in the channel with increase in gate bias. Eb(classical) is the bottom of the triangular well and the resonant 
level is shown by the horizontal line. 



We will now address the value of total transimssion in a ballistic MOSFET. The transmission is related to drain 
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current by (eqn. (II 5 9 [) ) . 



Id 




(209) 



where T is the total transmission from source to drain, fs and fo are the Fermi factors in the source and drain, and 
the factor of 2 accounts for spin. The main factors that determine the transmission probability are tunneling and 
scattering in the two dimensional potential profile, as an electron transits from source to drain. Traditionally simple 
theories of ballistic nanotransistors have assumed that the transmission from the source to drain to be integers in 
the above expression for current. The quantum mechanically computed transmission versus energy is shown in Fig. 
1171 The transmission increases in a step-wise manner, with the integer values at plateaus equal to the number of 
conducting modes in the channel. The steps turn on at an energy determined by an effective 'subband dependent' 
source injection barrier (E r i). That is, the maximum subband energy between the source and drain due to quantization 
perpendicular to the gate plane (x-direction of Fig. I12|) . The total transmission assumes integer values at an energy 
slightly above the maximum in 2D density of states as shown in the inset of Fig. [T7J Further, the transmission steps 
develop over a 50 meV energy window. In the case of MIT25 the current is predominantly carried at energies where 
the transmission is not an integer. 



FIG. 17: Transmission (+) and density of states (solid) versus energy at a spatial location close to the source injection barrier, 
at V g — 0V and Vd = IV. The peaks in the density of states represent the resonant levels in the channel. Inset: The density 
of states at three different y-locations and total transmission (+). The points y = -7 and nm are to the left and right of the 
location where the source injection barrier is largest (close to y = -4 nm). 



VIII. APPLICATION TO NANOTRANSISTORS WITH ELECTRON-PHONON SCATTERING 

The channel, scattering and screening lengths become comparable in transistors with diminishing channel lengths. 
Carrier transport is however not fully ballistic. Realsitic nanodevice modeling will involve phase-breaking scattering 
such that transport is between the ballitic and diffusive limits. In this regime, carriers are not thermally relaxed in 
the drain-end of the transistor, in contrast to long channel devices. The reflection of hot carriers from the drain-end 
towards the source-end, both above and below the source injection barrier should be explicitly included in models to 
compute the drive current. It is in this intermediate regime that the NEGF method has a distinct advantage over 
solving Schrodinger and Poisson equations self-consistently. 

To illustrate modeling of clcctron-phonon scattering in nanotransistors, we consider a prototype dual gate MOSFET 
(DGMOSFET) with a channel length of 25 nm and channel thickness (t c h) of 1.5 nm. The quantized energy levels in 
the channel due to quantization perpendicular to the gate (x-direction of Fig. I18[) are, 
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FIG. 18: Schematic of a Dual Gate MOSFET (DG MOSFET). Ex-s and Ex-d are the extension regions and the hatched region 
is the channel. The white region between the source / drain / channel and gate is the oxide. The device dimension normal to 
the page is infinite in extent. 



When t c h is small, the energy level separation is large and very few subbands are occupied in the highly doped 
extension regions. The lowest three quantized energy levels [eqn. (|210p ] are at 173, 691 and 891 meV above the 
bulk conduction band. The Fermi energy at the contact doping of 1E20 cm -3 is 60 meV above the bulk conduction 
band. As a result, electrons are injected only into one subband from the source-end at the operating voltage of 
Vd = V g = 0.6V. In this section, we discuss modeling of such MOSFETs using the mode space approach. The 
mode space approach consists of solving a one dimensional Schrodinger's equation for each subband. Inter subband 
scattering which can arise due to a change in electrostatic potential profile is neglected. Elcctron-phonon scattering 
between different subbands is included within the Born approximation. 

The three dimensional effective mass Hamiltonian is the same as eqn. (|193p . As the z direction is infinite, the wave 



function can be expanded as, &(x, y, z) 



(j)(x,y). Schrodinger's equation is then, 



E - 
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1 d 



4>{x,y) = 



(211) 



In the first step of the mode space approach, the eigen values [E n (y)\ at every cross section y along the source-drain 
direction is computed by solving for a subband dependent eigen value that varies with y, 

h 2 d 



2rwi dx 



1 d 

ink. dx 



V(x,y) 



y n {x,y) = E n (y)^ n (x,y) 



(212) 



n = {^,6}, where v is the quantum number due to quantization in the X-direction and 6=1,2 and 3 are the valley 
indices. 

In the second step, the Green's function equations for G, G n and G p are solved for each subband n: 



E 



h 2 k 2 



h 2 d 



l d 



2m" \ 2 dy \ m™ dy ' ™ ^ 



G n (y,y',k z ,E) 

-J dyi E n (y,y 1 ,k z ,E)G n (y 1 ,y',k z ,E) = S(y - y') , and (213) 



E 



h 2 k\ 

2m" 



h 2 d ( 1 d 
2 dy \ m™ dy 



G^(y,y',k z ,E) 

-J d Vl Z n (y 7yi ,k z ,E)G^ p) (yi,y',k z ,E) = J dy^ out Hy,y 1 ,k z ,E)Gi(y 1 ,y',k z ,m^ 

rriy and m z are the effective masses of silicon in the y and z directions that give rise to subband index n. Note that 
in Eqs. (|213[) and (|214[) . E n (y) is effectively an electrostatic potential for electrons in subband n. 

We designate the Green's/correlation functions G™ with a £ (empty), p, n and, respectively, self-energies S" with 
a G (empty), out, in. They can be written as, 



spa 



S" P , where 



n Ariel 



(215) 
(216) 
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E" c is the self-energy due to the leads. The phonon self-energy E" p consists of two terms, E" el due to elastic and 
due to inelastic scattering. Scattering between the lowest three subbands due to electron-phonon interaction 



is included and all other inter subband scattering is neglected. The self-energy due to leads is non zero only at the 
first (source) and last (drain) grid points. 

Assuming isotropic scattering, an equilibrium phonon bath and the self-consistent Born approximation, the self- 
energies due to electron-phonon scattering at grid point yi are jMah87l ] , 




G*,{ Vi ,E z ,E) 



(217) 
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l)GZ,( yi ,E z ,E + fiw v ) 



(218) 



and 
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inel.n 



(y h E) = y2Dy 




\n B (hu: n )G v n ,{y u E z ,E + hw n ) + {n B {hw n ) + l)G p n ,{y u E Z ,E- hu>r,)} 



(219) 



Here r\ represents the phonon modes, and the square of the matrix elements for phonon scattering are given by, 
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The contribution to elastic scattering is only from acoustic phonon scattering . The v alues of the deformation potential, 
Da, D gv and Df ri , and phonon frequencies ui grt and LUf v are taken from LunOOl ]. p is the mass density, k is the 
Boltzmann constant, T is the temperature and v is the velocity of sound, b and b' are indices representing the valley. 
The following scattering processes are included: acoustic phonon scattering in the elastic approximation and g-type 
intervalley scattering with phonon energies of 12, 19 and 62 meV. The imaginary part of the electron-phonon self- 
energy which is responsible for scattering induced broadening of energy levels and energetic redistribution of carriers 
are included but the real part is set to zero. 

In the numerical solution, N uniformly spaced grid points in the ^-direction with the grid spacing equal to Ay are 
considered. The discretized form of Eqs. (|2 13[) and (|214|) are then: 



A i ,iG n (y i ,y' i ,k z ,E) + A li+ iG n (y i+ i 1 y' ll k z ,E) + A i i - 1 G n (yi-i,y' i ,k z ,E) = , and 

Ay 



(222) 



Ai,iG%( yi , y' l ,k z ,E)+ Aij+tG^ivi+^yl, k z ,E) + A^G^y^^yl k z ,E) = 

^( yi ,E)GUy i ,y' i ,k z ,E) , (223) 



where, 
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The contact self-energies are: 



X n ,c(yi,k z ,E) = { —-^f gs{kz ,E) (226) 

^ n ,c{VN,k z ,E) = ( 2m f A ^ 2 ) 2 gd(fc z ,g) (227) 

H% c ( yi ,k„E) = -2lm(12 nt c(yi,k z ,E))f s (E) (228) 

^ c {y N ,k z ,E) = -2Im(£ n , c (y N ,k z ,E))f d (E) (229) 

Z°^( yi ,k z ,E) = -2Im(E n , c (yi,k z ,E))[l- f s (E)] (230) 

K*c(yN,k z ,E) = -21m(£ n , c (yN,k z ,E))[l-f d (E)] , (231) 



where y\ an y^ are the left (source-end) and right (drain-end) most grid points respectively, g s (k z ,E) and gd(k z ,E) 
are the surface Green's functions of the source and drain leads respectively, and f s and fd are the Fermi functions in 
the source and drain contacts respectively. 

The electron and current densities per energy given by Eqs. (|128[) and (|149[) can be simplified to 

n n ( yi ,k z ,E) = Gl{y uyi ,k z ,E) (232) 

J n (yi,k z ,E) = ~ 2m n A 2 l G n(Vi, Vi+i, k z ,E)- Gl{y i+U Vi , k z ,E)] . (233) 

n y " 

The total electron and current densities at grid point yi are given by, 

, x \frrW f dE f 1 s , 

= % J ~, J dE *7m nn{Vh E - E) (234) 

- 4 E ^/ f / ^.^^^^ , (235) 

where the prefactor of 4 in the above two equations account for two fold spin and valley degeneracies. The non 
equilibrium electron and current densities are calculated in both the channel and extension regions using the algorithm 
for G n presented in section [C"2l 

In solving the the Green's function and Poisson's equation, the applied bias corresponds to a difference in the 
Fermi levels used in the source and drain regions. The electrostatic potential at the left and right most grid points 
of the source and drain extension regions are calculated self consistently using the floating boundary condition that 
dV(y)/dy = 0. Poisson's equation is solved in two dimensions and the electron density is computed from Eqs. (12121) 
and (|232p using, 

n(xi,yi,k z ,E) = n n {y i ,k z ,E)\^ n {x i ,y i )\ 2 . (236) 

When does the mode-space approach fail? Ref. [Ven03l | has extensively analysed the regime of validity of the mode 
space approach. The mode space approach is valid when the wave function ^ n (x, y) at various y cross sections in eqn. 
(gUJ) satisfies ^ n{ d x y v) ~ 0. That is the shape of the wave function at each cross section should not change significantly 
along the transport direction, which means that intersubband scattering due to changes in potential profile is absent. 
This approximation seems to be valid until a channel thickness of 5 nm for silicon. 



A. Results 



Using the equations presented in the previous section, we show results illustrating the role of scattering along the 
channel length of a nanotransistor. First, we show that in devices where the scattering length is comparable to channel 
length, the nanotransistor drive current is affected by scattering at all points in the channel. Second, we show that 
when hot electrons enter the drain extension region of a nanotransistor, the drain extension region cannot be modeled 
as a series resistance. Instead, the drain extension should be included as part of the non equilibrium simulation region. 

To illustrate the role of scatter ing alo ng the channel, we calculate the drain current as a function of the right 
boundary of scattering (Yn-Scatt) Svi03| . Scattering is included only from the source end of the channel (-5 nm) to 



YR-scatt by setting the deformation potential in Eqs. (|220p and (|22ip to zero to the right of Yn_scatt- The scattering 
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FIG. 19: Drain current versus Ya-Scatt for two different scattering lengths. The channel length is 25 nm. 



lengths are decreased by a factor a by modifying the deformation potential in Eqs. (|220[) and (|22ip by an overall 
multiplicative factor of ^fa. 

The device considered has a channel length of 25 nm, body thickness of 1.5 nm, oxide thickness of 1.5 nm, doping 
of 1 E+20 cm -3 in the source and drain extension regions and an intrinsic channel. The scattering length due to 
clcctron-phonon interaction is 11 nm. We see from Fig. [TO] that scattering in the right half of the channel (0 to 12.5 
nm) is important and that the drive current degradation due to scattering in the right half of the channel is 30%. To 
understand the large reduction in drain current due to scattering in the right half of the channel, we plot the current 
density as a function of energy at various cross sections, J(Y, E). J(Y, E) shows the energetic redistribution of carriers 
along the channel. When the channel length is comparable to the scattering length, J(Y, E) in the right half of the 
channel is peaked in energy above the source injection barrier as shown in Fig. [20] (a). Scattering causes reflection 
of these energetic electrons toward the source. These reflected electrons lead to an increase in the channel electron 
density (classical MOSFET electrostatics). As the charge in the channel should be approximately C ox (Vg — Vg), the 
source injection barrier floats to a higher energ toy compensate for the refl ected electrons. The increase in source 
injection barrier and reflection leads to the decrease in drain current Lun97l | . 



To gain further insight into the role of carrier relaxation, we now discuss the case when the scattering length is 
2.2 nm, which is five times smaller than in the previous discussion. We see from Fig. [151 that scattering in the right 
half of the channel now decreases the drive current by a smaller amount of 15%, when compared to the case with 
Lscatt = 11 nm. As the channel length of 25 nm is much larger than 2.2 nm, multiple scattering events now lead to an 
energy distribution of current that is peaked well below the source injection barrier in the right half of the channel as 
shown in Fig. [^D] (b) . The first moment of energy with respect to the current distribution function, which is defined 
as the ratio of / dEEJ(Y, E) and total current is also shown in Fig. [20] When the scattering length is much smaller 
than the channel length, the carriers relax classically such that the first moment (J dEEJ(Y, E)) closely tracks the 
potential profile as seen in Fig. [20] (b) . Thermalized carriers that are reflected in the right half of the channel can no 
longer reach the source injection barrier due to the large barrier to the left, and so contribute less significantly to the 
charge density. Thus, explaining the diminished influence of scattering in the right half of the channel relative to the 
left half of the channel. 

To further demonstrate the use of NEGF simulations, we study the role of assuming that the extension regions can 
be modeled as a classical series resistance. Within the classical serie s resista nce picture the current with scattering 
(Ip att ) can be related to the current without scattering (/™ oscatt ) by [Tau98l ] . 

I s D catt (.V D ) ~ I~ U (V D - SV D ) , (237) 

where we have assumed that the source extension region and device do not experience scattering. The potential drop 
in the drain region within the classical series resistance picture is 8Vd — Ifj att (VD)RD- In Fig. [2T] the values of 
the drain current in the ballistic limit, using the classical series resistance picture and with the NEGF method are 
marked. Clearly, the series resistance picture underestimates the deterimental nature of scattering in the drain end. 
The physics of the large reduction in drain current was discussed in the context of Fig. [20J When scattering in 
the channel does not effectively thermalize carriers, the current distribution is peaked at energies above the source 



35 




Y (nm) Y (nm) 



FIG. 20: Solid lines represent J(Y,E) for Y equal to -17.5, -12.5 -7.5, -2.5, 2.5, 7.5, 12.5, 17.5 nm respectively, when scattering 
is included every where in the channel. The dashed lines are the first resonant level (Ei) along the channel. The dotted lines 
represent the first moment of energy with respect to the current distribution function, J dEEJ(Y, E). (a) and (b) correspond 
to L SC att = 11 and 2.2 nm respectively. 
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FIG. 21: Ballistic Id{Vd) with the drive currents obtained in the ballistic limit, with the series resistance picture and NEGF 
calculations marked. 

injection barrier, upon carriers exiting the channel. Scattering in the drain extension region then causes reflection of 
electrons toward the source-end. As a result, the source injection increases so as to keep the electron density in the 
channel approximately at C ox (Vg — Vs). The drain current decreases dramatically as a result of the increase in source 
injection barrier height. 

IX. DISCUSSION AND SUMMARY 

Our objectives in this chapter has been to: (i) review the underlying assumptions of the traditional, semiclasssical 
treatment of carrier transport in semiconductor devices, (ii) describe how the semiclassical approach can be applied 
to ballistic transport, (iii) discuss the Landaucr-Buttikcr approach to quantum transport in the phase coherent limit, 
(iv) introduce important elements of the non-equilibrium Green's function approach using Schrodingcr's equation as 
a starting point and (v) finally demonstrate the application of the NEGF method to the MOSFET in the ballistic 
limit and with electron-phonon scattering. 

It is appropriate to make a few comments about the computational burden of the various transport models. One 
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reason that device engineers continue to use drift-diffusion simulations rather that the more rigorous Monte Carlo 
simulations is the enormous difference in computational burden. For semiclassical transport, the fundamental quantity 
is the carrier distribution function, f(f,k,t). To find f(f,k), we must solve the BTE, which is a six-dimensional 
equation. The difficulty of solving this six-dimensional equation is one reason that engineers continue to rely on 
simplified models. For quantum transport, we can take the Green's function G n (r, f", E) as the fundamental quantity. 
The Green's function is a correlation function that describes the phase relationship between the wavefunction at r and 
f 1 for an electron injected at energy E. The quantum transport problem is seven dimensional, which makes it much 
harder than the semiclassical problem. We can think of r and f* as analogous to r and k in the semiclassical approach, 
but there is no E in the semiclassical approach. The reason is that for a bulk semiconductor or in a device in which the 
potential changes slowly, there is a relation between E and k, as determined by the semiconductor bandstructure E(k). 
When the potential varies rapidly, however, there is no E(k), and energy becomes a separate dimension. Analysis 
of electronic devices by quantum simulation is however becoming practical because device dimensions are shrinking, 
which reduces the size of the problem. Quantum simulations are also essential to accurately model devices whose 
dimensions are comparable to the phase breaking length, and rely on tunneling and wave interference for operation. 
The resonant tunneling diode is the most successful example in this category. 
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APPENDIX A: DERIVATION OF EQN. ([60]) 

Eq. can be expanded as, 

G LD = -A'-lA' LD G DD (Al) 
Grd = —A'£ r A' rd Gdd (A2) 
A'dlGld + A' dd Gdd + A dr Grd = I ■ (A3) 

Substituting Eqs. (|A"Tj) and (|X2|) in cqn. (|X3|) . we have, 

[A'dd ~ a dl a ll a 'ld ~ A' dr A' rr ~A' rd \Gdd = I ■ (A4) 

Noting the sparsity of A' LD and A' RD , it follows that only the surface Green's functions of the Left and Right 
leads, 

= a 'll 1a and 9R n+1 , n+1 = A' RR - 1A (A5) 
are required in eqn. (|A4[) . Then, we can rewrite eqn. (|A4[) in terms of lead self-energies as, 

AG DD = [EI -H D - ^i ead ]G DD = I, (A6) 
where S/ ea d has been defined in eqn. (f52"| and (f53"| . The reader can verify that gr J _ 1 _ 1 = e +lA;,a i ; ~ 1 and gn n+1 n+1 = 

g+ifcrO^— 1 

APPENDIX B: DYSON'S EQUATION FOR LAYERED STRUCTURES 
(Note: Matrix A of this section is equivalent to Matrix Ajyjy of section IVl]) 

Partition the device layers into two regions Z and Z' as shown in Fig. [221 Dyson's equation is a very useful method 
that relates the Green's function of the full system Z + Z' in terms of the subsystems Z' and the coupling between 
Z and Z 1 . We will see below that from a computational point of view, Dyson's equation provides us with a systematic 
framework to calculate the diagonal blocks of G and G n without full inversion of the A matrix. The re ader sho uld 
note that Dyson's equation has a significantly more general validity than implied in our application here Mah9dl | . 
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FIG. 22: Scheme of device for application of Dyson's equation by splitting the device in two parts. 



1. Dyson's equation for G 



The Green's function equation over the device layers [eqn. (11221) ]. 

AG = I 



can be written as, 



The solution of cqn. (jBll s. 



where, 



Az,z Az t z' \ I Gz,z Gz,z' 
Az'z Az'.z' ) \ Gz'.z Gz'Z' 



G = G° + G°UG 
= G° + GUG° 



I O 
O I 



G = 



Gz,z Gz,z' 
Gz'Z Gz'z 1 



qO _ ( G z,z O 



O A 



and U = 



Z'Z> 



The Hcrmitean conjugate Green's function (G>) is by definition related to G by 

G t = G^ + G^U^G" 1 
= G t0 + G^U^G^ . 

Eq. (|B3|) is Dyson's equation for the Green's function. 



O -A z ,z> 
Az'z O 



(Bl) 



(B2) 



(B3) 
(B4) 
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(B6) 
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2. Dyson's equation for G n 

The electron correlation function equation over the device layers [eqn. (|131| ]. 

AG n = S™G f (B8) 



can be written as, 



The solution of cqn. (|B8[) is, 



Az.z Az,z> j | G z z G z z , j _ | ^z)z' j ( G z.z G z .z> j j-gg-j 

A Z ',z A z >,z' J \ G Z 'z G z',z< j V ^z\z ^z'.z' J \ G z>,z G z>,z> J 



G n = G°UG n + G°Y, in G\ (BIO) 
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where G° and U have been defined in Eqs. (|B5[) . Functions G n and G^ are readily defined by eqns. (|B8[) and (|B7|) . 
respectively. Using G 1 * = G t0 + G^WG\ eqn. (|B10p can be written as 

G" = G"° + G™°f/ t G t + G°f/G™ (Bll) 

= G n0 + GUG ,M + G n U i <G f0 , (B12) 

where G"° = G°S m G t0 . (B13) 



APPENDIX C: ALGORITHM TO CALCULATE G AND G 



(Note: Matrix A of this section is equivalent to Matrix Add of section IVlf 



Why algorithm: A typical simulation of a nanoelectronic device consists of solving Poisson's equation self- 
consistcntly with the Green's function equations. The input to Poisson's equation is the charge density which is 
obtained via integrating over energies the elements, G™ q (E), from the main diagonal of the electron correlation func- 
tion. The index q here runs over the layers of the device. In order to calculate the current density one requires the 
elements, G q q+1 (E), from the diagonals adjacent to the main diagonal. That is, we do not require the entire G™ 
matrix in most situations. Same goes for the hole correlation function G p and the Green's function G. 

Provided that N x is the dimension of the Hamiltonian of each layer and N y is the total number of layers, the size 
of the matrix A equals N x N y The operation count for the full matrix inversion G = A^ 1 is proportional to (N x N y ) 3 . 
The computational cost of obtaining the diagonal elements of the G™ matrix at each energy is approximately N^Ny 
operations if G n = GE m G^ is used. Therefore it is highly desirable to find less expensive algorithms that avoid full 
inversion of matrix A and take advantage of the fact the diagonal elements of Green's functions. Another reason to 
prefer such algorithms is the memory storage. If one had had to retain the whole matrix G in the memory, it might 
had required using the RAM or the hard drive instead of on-chip cache. That would have significantly slowed down 
the calculations. 

One such algorithm which is valid for the block tridiagonal form of matrix A is presented in this section. The 
operation count of this algorithm scales approximately as N^N y . The dependence on N% arises because matrices of 
dimension of the sub Hamiltonian of layers should be inverted, and the dependence on N y corresponds to one such 
inversion for each of the N y layers. 

The algorithm consists of two steps. In the first step, the diagonal blocks of the left connected and full Green's 
function are evaluated (section [C In the second step, these results are used to evaluate the diagonal blocks of the 
less-than Green's function (section [G 2[) . 



1. Recursive algorithm for G 



(i) Left-connected Green's function (Fig. |23|) : 

The left-connected (superscript L) Green's function g Lq is defined by the first q blocks of eqn. (|122p by 

g Lq = /,,,. (CI) 

where we introduce a shorthand I q ^ q = 2i :g ,i :g The matrix g Lq+1 is defined in a manner identical to g Lq except that 
the left-connected system is comprised of the first q + 1 blocks of eqn. (|122[) . In terms of eqn. (|B2[) . the equation 
governing g Lq+1 can be expressed via the solution g Lq by setting Z = 1 : q and Z' = q + 1. Using Dyson's equation 
[eqn. (|P33|) ] . we obtain 

9g+tg+l = (Aq+l,q+l ~ A q+ l, q g q *A q , q+ l) . (C2) 

Note that the last element of this progression g LN is equal to the fully connected Green's function G, which is the 
solution to eqn. (|122|) . 

(ii) Full Green's function in terms of the left-connected Green's function : 

Consider the special case of eqn. (|B2[) in which Az,z = Ai-.q,\:qi Az\z' = Aq+i:N,q+i:N and Az,z' = ^4i: g , g +i:iv- 
Noting that the only nonzero element of Ai- qtq +x;N is A q ^ q +\ and using eqn. (|B3[) . we obtain 

Gq.q = 9 q q q + 9q, q q (^-q,q+lGq+l,q+l^-q+l,q) 9q, q q (C3) 
= 9q, q q - 9q q q Aq,q+lG q+ l. q . (C4) 
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terms of e 1 *, p rti * using Dyson's equation. 



FIG. 23: Illustration for the relation between the left-connected Green's functions for adjacent layers. 



The equations for the adjacent diagonals are obtained similarly 

Gq+i.q = -G q+ i^ q+ iA q+ i^ q g^ q q , (C5) 
G q , q +\ = -g q q q A q ^ q+ iG q+ i }q+ i . (C6) 

Both G q ^ q and G g +i, g are used in the algorithm for electron density, and so storing both sets of matrices is necessary. 
Making use of the above equations, the algorithm to obtain the three diagonals of G is 

i n Ll _ A-l 

2. For q=l,2, ...,7V - 1, compute eqn. (jC2|) . 

3. For q= 1,2, ...,7V, compute (g^f . 

4. Gn.n = g q q q - 

5. For q = TV - 1,7V - 2, 1, compute eqns. ([C5]t . (fC6|) and ([C4l (in this order). 

6. For g = l,2, ...,7V, compute (G q . q+ \)^ and {G q+ \^ . 



2. Recursive algorithm for G 



(i) Left-connected G n (Fig. 

The function g nLq is the counterpart of g Lq , and is defined by the first q blocks of eqn. (|13ip : 

A 1:g>1: , = E*? gjl: , . (C7) 

griLq+i - 1S ^ e fi nec j m a manner identical to <7™ L|? except that the left-connected system is comprised of the first q + 1 
blocks of eqn. Q13ip . 

The equation governing g nL i+ x follows from eqn. (|B9[) by setting Z = 1 : q and = q + 1. Using the Dyson's 
equations for G [eqn. (|B3j) ] and G n [eqn. ([Blip ], is recursively obtained as 



nLq+l _ Lq+l Ty,i n . i n ] Lg+lf fpo^ 

where cr* 1 ^ = A q +i. q g^ q A q +1 . Eqn. (|C8j) has the physical meaning that ff^+i^+i has contributions due to an 
effective self-energy due to the left-connected structure that ends at q, which is represented by <J q \ 1 and the diagonal 
self-energy component at grid point q + 1 of eqn. (|13ip ). 
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(ii) Full electron correlation function in terms of lcft-conncctcd Green's function : Consider eqn. (|B9[) such that Az = 
Ai:q,i :q , A' z = A g +i : jv,g+i:iV and Az.z> = ^-l-.q.q+i-.N ■ Noting that the only nonzero element of A 1:q ^ q+1 . N is A q>q+ i 
and using eqn. (|Bllj) . we obtain 

rm n nL 1 — n nL 1 A^ — n^l A C n /pn\ 

q,q ~~ ilq.q iJq,q ■ fl q,q+l lJ q+l,q yq,q- f± q,q+l'~ r q+l.q • 

Using eqn. (|B12|) . G q+1 q can be written in terms of G q+1 q+1 and other known Green's functions as 

Gq+i,q = ^G q+ i^q + iAq + i iq g q L q q — G q+lq+1 A q+1 q g\ L q ■ (CIO) 
Substituting eqn. (|C10|) in eqn. (|C9[) and using Eqs. (|B3|1 and (|B4|) . we obtain 

rm n n Lq i „Lq ( a rm /it 1 „\Lq _ „ n Lq a] rA \ n A r, nL q 

~ ilq.q ilq,q \ J± q,q+l'~ r q+l.q+l I1 q+l,q J !Jq,q ilq,q -™-q,q+l ^q+l.q > ^9,9+1 ■ e± q+l,q!fq,q 
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(Cll) 



The terms inside the square brackets of eqn. (|C11[) are Hermitian conjugates of each other. In view of the above 
equations, the algorithm to compute the diagonal blocks of G n and G p is given by the following steps: 

1 n nLl _ Llyin ilt 

2. For q=l,2, ...,7V - 1, compute eqn. l(C8|) . 

o rm ^nLN 

°- U N,N — 9NN ■ 

4. For q = N - 1,N - 2, 1, compute eqns. (|Cllj) and (jClOjl . 

5. Use G" g+1 = (G™ +lg ) . 

6. Use G» = i(G-Gt) - G" . 

The above algorithm is illustrated by a Matlab code in Section |Pl 
Challenging problem: The algorithm presented here solves for the three block diagonals of G, G", and G p . Each 
of n blocks on the main diagonal corresponding to a layer of the device. All blocks in the three diagonals is treated 
as a full matrix. It is highly desirable to find a more efficient algorithm that finds only the diagonal elements within 
each block rather than complete blocks. 

APPENDIX D: CODE OF THE RECURSIVE ALGORITHM 

The listing of a Matlab code recursealg3d.m is provided here for illustration of the algorithm described in Section 



f unct ion [Grl , Grd , Gru , Gnl , Gnd , Gnu , Gpl , Gpd , Gpu , grL , ginL] = recur sealg3d (Np , Al , Ad , Au , Sigin , Sigout ) 

°/ function [Grl, Grd, Gru, Gnl, Gnd, Gnu, Gpl, Gpd, Gpu] = recursealg3d(Np , Al , Ad, Au, Sigin, Sigout) 

°/ recursive algorithm to solve for the diagonal elements of 

°/« the Non-equilibrium Green's function 

'/, "1" = lower diagonal = [G(2,l); ...; G(end,end-1) ; 0] 

% "d" = main diagonal = [G(l,l); ...; G(end,end)] 

7, "u" = upper diagonal = [0; G(l,2); G(end-l.end)] 

°/« Grl, Grd, Gru = retarded Green's function 

°/ Gnl, Gnd, Gnu = electron Green's function 

°/ Gpl, Gpd, Gpu = hole Green's function 

7, Np = size of the matrices 

°L Al,Ad,Au = matrix of coefficients 

°/ Sigin = matrix of in-scattering self-energies (diagonal) 

% Sigout = matrix of out-scattering self -energies (diagonal) 

°/ Dmitri Nikonov, Intel Corp. and Siyu Koswatta, Purdue University, 2004 

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 

f lag_Gp = 'no ' ; 
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Al_cr = conj (Au) 
Ad_cr = conj (Ad) 
Au_cr = conj (Al) 
grL = zeros (l,Np); 
ginL = zeros (l,Np); 
gipL = zeros (l,Np); 
Grl = zeros (1, Np-1) ; 
Grd = zeros (l,Np) ; 
Gru = zeros(l ,Np-l) ; 
Gnl = zeros(l ,Np-l) ; 
Gnd = zeros (1 ,Np) ; 
Gnu = zeros(l ,Np-l) ; 
Gpl = zeros(l ,Np-l) ; 
Gpd = zeros(l,Np); 
Gpu = zeros(l ,Np-l) ; 
grL(l)=l/Ad(l) ; 
for q=2:Np 



7 Hermitean conjugate of the coefficient matrix 

7 initialize left-connected function 

7 initialize left-connected in-scattering function 

7 initialize left-connected out-scattering function 

7 initialize the Green's function 



7 initialize the electron coherence function 



7 initialize the hole coherence function 



Grd(Np)=grL(Np) 
for q=(Np-l) :-l 



7 step 1 

7 obtain the left-connected function 
grL(q)=l/(Ad(q)-Al(q-l)*grL(q-l)*Au(q-l)) ; 

end 

gaL = conj (grL) ; 7 advanced left-connected function 

7 step 2 

1 

Grl(q)=-Grd(q+l)*Al(q)*grL(q) ; 7 obtain the sub-diagonal of the Green's function 

Gru(q)=-grL(q)*Au(q)*Grd(q+l) ; 7 obtain the super-diagonal of the Green's function 

Grd(q)=grL(q)-grL(q)*Au(q)*Grl(q) ; 7 obtain the diagonal of the Green's function 

end 

Gal = conj (Gru) ; 

Gad = conj (Grd) ; 7 advanced Green's function 

Gau = conj (Gal) ; 

ginL(l)=grL(l)*Sigin(l)*gaL(l) ; 7 step 3 

for q=2:Np 

sla2 = Al(q-l)*ginL(q-l)*Au_cr(q-l) ; 

prom = Sigin(q) + sla2; 

ginL(q)=grL(q)*prom*gaL(q) ; 7 left-connected in-scattering function 

end 

Gnd(Np)=ginL(Np) ; 7 step 4 

Gnd = real (Gnd) ; 
for q=(Np-l) :-l:l 

Gnl(q) = - Grd(q+l)*Al(q)*ginL(q) - Gnd(q+l)*Al_cr(q)*gaL(q) ; 

7 obtain the lower diagonal of the electron Green's function 

Gnd(q) = ginL(q) + grL(q) *Au(q) *Gnd(q+l) *Al_cr (q) *gaL(q) ... 

- ( ginL(q)*Au_cr(q)*Gal(q) + Gru(q) *A1 (q) *ginL(q) ); 

end 

Gnu = conj (Gnl) ; 7 upper diagonal of the electron function 

switch flag_Gp 
case 'yes' 

gipL(l)=grL(l)*Sigout(l)*gaL(l) ; 7 step 3 

for q=2:Np 

sla2 = Al(q-l)*gipL(q-l)*Au_cr(q-l) ; 

prom = Sigout(q) + sla2; 

gipL(q)=grL(q) *prom*gaL(q) ; 7 left-connected in-scattering function 

end 

Gpd(Np)=gipL(Np) ; 7 step 4 

Gpd = real (Gpd) ; 
for q=(Np-l) :-l:l 

Gpl(q) = - Grd(q+l)*Al(q)*gipL(q) - Gnd(q+1) *Al_cr (q) *gaL (q) ; 

7 obtain the lower diagonal of the hole Green's function 

Gpd(q) = gipL(q) + grL(q) *Au(q) *Gpd(q+l) *Al_cr(q) *gaL(q) ... 
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- ( gipL(q)*Au_cr(q)*Gal(q) + Gru(q)*Al(q)*gipL(q) ); 

end 

Gpu = conj (Gpl) ; °L upper diagonal of the hole function 

case 'no' 
Gpl = i*(Grl-Gal) - Gnl; 

Gpd = real (i* (Grd-Gad) - Gnd) ; % hole Green's function 

Gpu = i*(Gru-Gau) - Gnu; 

end 

jnzer = f ind(Gnd<0) ; 
Gnd(jnzer) = 0; 
jpzer = f ind(Gpd<0) ; 
Gpd(jpzer) = 0; 
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